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INTRODUCTION 


Aerosols  in  the  atmosphere  emit  thermal  radiation  with  intensities 
which  are  functions  of  their  temperature  and  wavelength.  The  thermal 
emission  by  airborne  aerosols  is  often  a significant  factor  in  deter- 
mining the  intensity  of  infrared  radiation  transmitted  through  the 
atmosphere.  For  example,  in  determining  vertical  temperature  profiles 
from  radiance  measurements  made  from  satellites,  the  thermal  emission 
from  clouds  within  the  field  of  view  of  the  instrument  must  be  con- 
sidered. Another  problem  of  considerable  interest  to  the  Army  in 
which  aerosol  thermal  emission  may  be  important  is  in  battlefield 
scenarios  in  which  smoke  clouds  may  interfere  with  weapon  systems 
dependent  upon  infrared  viewing  devices. 

The  approach  used  here  in  calculating  the  thermal  radiation  is  .analo- 
gous to  the  "doubling"  technique  used  in  computing  aerosol  multiple 
scattering  effects.  A summary  of  the  doubling  technique  is  presented, 
followed  by  a discussion  of  its  extension  to  thermal  emission  calcu- 
lations. A verification  is  discussed  in  which  results  are  compared 
with  some  published  results  of  thermal  emission  from  clouds.  Also 
presented  are  the  results  of  theoretical  studies  determining  the 
effects  of  surface  radiation  passing  through  a simulated  cloud  of 
smoke.  The  computer  program  written  for  the  Univac  1108  used  in 
these  studies  is  documented  in  the  appendix.  This  program  is  capable 
of  computing  water  vapor  continuum  absorption  in  the  8-  to  12- 
micrometer  region  using  a semiempirical  technique  [1], 

SUMMARY  OF  THE  DOUBLING  TECHNIQUE 

Consider  a horizontally  infinite  layer  of  homogeneous  aerosol  material 

_T 

whose  transmittance  in  the  absence  of  multiple  scattering  is  e , 
where  t is  the  optical  depth  of  the  layer.  Radiation,  in  the  form  of 

9 

a collimated  beam  of  irradiance  ttF  (watts/meter  ),  is  assumed  to  be 
incident  from  above  upon  the  top  of  the  layer  at  angle  0q  and  q 

(see  figure  1).  It  is  conventional  to  represent  6q  by  its  cosine, 

Wq  * cos  0q,  so  that  directions  are  specified  by  p and  d>  rather  than 

by  0 and  $ . To  avoid  using  negative  values  of  p , a convention  is 
adopted  in  which  "downward"  going  radiation  is  specified  by  ( — p , <{> ) 
and  "upward"  going  radiation  by  (p,4>);  thus,  the  p's  are  always  re- 
garded as  positive  numbers.  With  these  conventions,  the  incident 
beam  direction  is  specified  as  (-p^,^). 

2 

Now,  let  l(0,p,<f>)  be  the  radiance  (watts/meter  -steradian)  reflected 
in  the  direction  (p,<f>)  from  the  top  of  the  layer.  In  addition,  let 
the  radiance  transmitted  out  from  the  bottom  of  the  layer  (at  optical 
depth  t)  be  denoted  by  I(t,-p,<|>).  In  the  following,  the  subscript  v 
is  suppressed,  although  it  should  be  understood  unless  explicitly 
stated  otherwise. 


3 


(p  negative) 


g.  1.  Angular  definitions  for  multiple  scattering  calculations 
(above)  and  illustration  of  radiance  vectors  in  the  planes 
$ -41  =0®  and  ♦ = 180®. 


The  doubling  approach  to  radiative  transfer  problems  seeks  the  calcu- 
lation of  two  "operators"  S and  T with  the  properties: 


I (0 ; y , <t>)  “ “J  S(t;u,<I>; > d>0>  * <ttF)  , (1) 

I (t y > 4>)  = T(T;-y,(();uQ,0Q)  • (ttF) 


(nF)  -Tl/P0 

+ -^e  6(p0  - u)«4  - 4»0)-  (2) 

The  second  term  of  equation  (2)  represents  the  attenuated  direct  beam 
while  the  first  term  represents  diffusely  transmitted  radiance.  In 
subsequent  expressions,  the  factor  ttF  is  taken  to  be  unity. 

The  quantities  S/  (4irp)  and  T/ (4up ) can  be  regarded  as  operators  whose 

operands  are  the  incident  irradiance.  (Alternatively  S/P  and  T/P 
operate  on  incident  radiance.)  It  may  be  convenient  at  times  to  sup- 
press the  u ^ factors  by  writing  S for  S/p  and  T for  T/p.  The  meanings 
of  S and  T can  be  seen  as  follows:  let  1^  = I(0;p,4>)  , the  upward  in- 
tensity at  t = 0;  Iq  = I(0;-Pq,<J)q)  , the  downward  intensity  at  x = 0, 
and  1^  = I(T^;-P»<I>) » the  downward  intensity  at  t = x . With  these 
definitions,  one  can  write 

i;  - S ij  and  ■ T ij, 

Multiple-scattering  calculations  using  the  doubling  approach  are  based 
upon  assumptions  that:  (1)  one  can  find  expressions  for  S and  T which 
are  valid  at  very  small  optical  depth,  that  is,  under  single-scattering 
conditions,  and  (2)  one  '••'n  define  a procedure  for  finding  the  overall 
scattering  and  transmit  functions  for  a combination  of  two  layers 
whose  individual  S and  inctions  are  already  known.  To  this  end, 
the  layer  is  broken  up  i..to  a number  of  layers  such  that  eventual 
summation  of  those  layers  will  have  the  same  optical  depth  as  the 
original  layer. 

To  determine  initial  values,  one  starts  with  a layer  of  very  small 

—20 

optical  thickness  (~2  ) such  that  the  initial  S and  T values  are 

essentially  given  by  the  phase  function  [2,3].  Once  the  initial 
values  have  been  determined,  the  values  of  S and  T for  a layer  of 
twice  that  starting  thickness  (optical  depth)  are  found  by  "stacking" 
two  such  layers  together.  The  stacking  process  is  then  repeated  until 
the  desired  final  x value  is  reached.  This  general  procedure  is  widely 
known  as  the  "doubling"  method  [3,4], 


5 


It  is  assumed  that  S and  T can  be  expanded  in  a Fourier  series  over  $ 
and  (}>q  as  follows :+ 

N , ) r 1 

S(T;u,4>;u0.<f>0)  *=  m~0  S m;(T;y,P0)  cos  ^m(<l>  - <)>0)J  , (3) 

and 

T(t ;y » i u q » ^ 0>  (t  ;vj,ViQ)  cos 

Furthermore,  since  the  calculations  are  numerical  (u  and  are  calcu- 
lated for  a discrete  number  of  angles  and  ik),  the  operators  S and 
T may  be  expressed  as  (operator)  matrices  with  the  following  notational 
changes: 

S(m)(T,li,P0)  + S®  (t)  - S"  , (5) 


T(m)(Ttli,y0)  - (t>  - Tj  . (6) 

Initial  or  starting  values  of  and  are  found  from  the  following 
expressions  [4]: 

s>o>  “ (ai  + aj)_1  {1_exp  [”(ai  + °j)To] } P!j  . (7) 

Tij(xo)  = (ai  + aJ)"1  {exp("°jTo)  ' exp(-aiTo)}  pij  . (8) 

where 


^For  problems  of  azimuthal  symmetry,  such  as  the  emissions  of  a black- 
body,  only  the  term  m ■ 0 Is  required. 


= 2 - « 

(12) 

0,m, 

({,  - m)  ! ~ 

= — ; — 0). 

(13) 

(H  + m)  ! *• , 

and  Pm  means  P?(M.),  an  associated  Legendre  polynomial  evaluated  at 

Hi  *•  1 

p = p . t is  the  optical  depth  of  the  initial  layer  and  is  taken 
i 0 _o  o 

typically  to  be  2 

For  completeness,  the  are  determined  in  the  single  scattering 
routines  (by  Gauss-Legendre  quadrature)  from: 


2 a 


TT 


p(cos  0)Pri  (cos  8)sin  0 d0 


(14) 


in  which  p(cos  0)  is  the  "phase  function"  and  0 is  the  scatter- 
ing angle: 


Once  and  T™.  have  been  calculated  for  an  initial  layer,  one  proceeds 

to  combine  two  such  layers  by  some  "process."  The  algorithm  used  here 
for  doing  that  will  be  reviewed  next.  For  simplicity,  again,  all  sub- 
scripts and  functional  dependencies  will  be  suppressed:  S will  rep- 
resent S^(t),  1q  will  represent  I(0;-Pq,<J>  ),  etc. 


Consider  two  layers  identified  by  the  numerals  1 and  2,  with  1 being 
the  "upper  layer."  The  upper  layer  is  to  be  characterized  by  S and 

* * ^ 

T for  illumination  from  above,  and  S , T..  for  illumination  from  below. 
1 1 1 * 


In  a similar  fashion  the  lower  layer  is  characterized  by  S2>  T2,  S2, 
T*  (asymae tries  for  illumination  from  above  or  below  are  "permitted") 


Let 

it  be  the  downward  intensity  incident  on  the  top  of  the  upper 
layer, 

1^  be  the  upward  intensity  reflected  from  the  top  of  the  upper 
layer, 

and 

I2  be  the  downward  intensity  at  the  bottom  of  the  lower  layer. 
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Reference  to  Figure  2 may  be  helpful. 

Using  these  definitions,  Barkstrom  and  Querfeld  [4]  have  shown  that 

I*  = SI*  and  I*  = TI* , 

0 0 2 0 

where 

T = T2(l  - S^)’1^  (16) 

and 

S = S1  + T*S2(1  - S*S2)"1T1  . (17) 

The  above  expressions  S and  T were  written  in  such  a way  as  to  permit 
stacking  nonidentical  layers,  if  desired.  Note  that  S^,  T^,  etc.,  are 

matrices  and  that  equations  such  as  (17)  imply  matrix  multiplication 
as  well  as  the  Fourier  summations  indicated  in  equations  (3)  and  (4). 

INCLUSION  OF  THERMAL  EMISSION 

For  a very  thin  layer  at  a temperature  t in  which  single  scattering 
dominates,  the  emitted  intensity  can  be  well  approximated  by  the 
expression 

I(^  = (1  - <0(1  - e‘T0/y)B  (t),  (18) 

vE  0 v 

where  is  the  optical  depth  of  the  layer,  o>q  is  the  single  scattering 
albedo,  and  Bv(t)  is  the  Planck  function  at  frequency  v.  If  one  re- 
gards the  intensity  emitted  by  one  layer  as  part  of  that  incident  upon 
another  one,  it  may  be  possible  to  assemble  an  emission  operator 
analogous  to  S or  T.  A method  for  doing  that  will  be  described  next. 

Recall  now  the  two  layers  utilized  in  the  previous  discussion.  Ignoring, 
for  the  present,  any  external  sources  of  radiation,  one  can  introduce 
some  additional  definitions  and  quantities.  Define 

I*  as  the  upward  directed  thermal  intensity  found  at  the 
top  of  the  upper  layer  in  the  absence  of  the  lower 
layer, 

1^2  as  the  downward  directed  thermal  emission  intensity 
at  the  bottom  of  the  upper  layer  in  the  absence  of 
the  lower  one. 
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Ik.- 


f 


and 


2E 


I1 

2E 


3 


3 


-i 


M 


as  the  upward  self-intensity  at  the  top  of  the  lower 
layer  in  the  absence  of  the  upper  one, 

as  the  downward  directed  self-emission  intensity  at 
the  bottom  of  the  lower  layer  in  the  absence  of  the 
upper  layer, 

as  the  total  upward  thermal  emission  intensity  at  the 
top  of  the  upper  layer  when  the  lower  layer  is  also 
present, 

as  the  total  downward  intensity  seen  from  the  bottom 
of  the  lower  layer  when  the  upper  layer  is  also  present, 

as  the  total  upward  directed  intensity  in  a hypo- 
thetical interstitial  space  between  the  two  layers 
when  both  exist, 

as  the  total  downward  intensity  in  the  interstitial 
space. 


The  quantity  1+  is  the  sum  of:  (1)  the  intensity  emitted  directly  by 

L 

the  upper  layer,  and  (2)  that  portion  of  1^  transmitted  from  below  by 
the' upper  layer  itself.  Thus,  one  can  immediately  write  the  expressions: 


1+  = 1+  + T 1+  and  1 + 

E IE  1 M E 


I'*'  + T I+  . 
2E  2 M 


(19*,b) 


The  quantity  I*  is  the  sum  of  I*  and  that  portion  of  I*  which  is  re- 
M 2E  M 

fleeted  by  the  lower  layer,  and  I+  is  similarly  related  to  I*  and 

M IE 

I+.  Therefore, 

M 


*4 


XJe  + S2XM and  XA  - XJe  + SA 


(20a, b) 

The  last  two  equations  can  be  combined  to  give,  for  example, 

:M  * (1  - SIS2>'1[IlE  + Sl 

and  (21a, b) 

■ a - S!V1[12E  + S2  11E> 


which  can,  in  turn,  be  substituted  into  equation  (19b)  to  yield: 

* 


1 + ■ iL  + (T0QS, ) ll_  + (T.Q)I^ 
E 2E  2H  1 2E  ' 2X  IE 

10 


(22) 


where  Q stands  for  (1  - S^S^)  . In  a similar  way  one  finds  that: 

li  ' *k  + (IJS2Q)I1E  + T1Q  *k  ' (23> 

Equations  (22)  and  (23)  provide  for  combination  of  the  self-emission 
intensities  of  two  layers  whose  individual  properties  are  known  but 
not  necessarily  identical.  Thus,  they  may  be  used  to  combine  the 
thermal  emissions  from  two  layers  which  differ  from  one  another  in 
temperature  as  well  as  in  passive  optical  properties.  Equation  (18) 
may  then  be  used  to  provide  starting  values  for  successive  applications 
of  (22)  and  (23). 

The  formal  procedures  derived  above  may,  in  principle,  be  useful  for 
the  construction  of  thermal  emission  intensities  for  a cloud  composed 
of  nonidentical  layers  of  arbitrarily  small  thickness  and  continuously 
varying  temperature.  That  sort  of  procedure  is  not  practical,  however, 
because  of  the  inordinately  large  number  of  stackings  which  would  be 
required  to  attain  a total  thickness  of  real  interest.  To  apply  equa- 
tions (22)  and  (23)  to  a more  realistic  case  (insofar  as  numerical 
computations  are  concerned) , a cloud  must  be  divided  into  a small 
number  of  layers  (~  20)  within  each  of  which  temperature,  droplet 
size  distributions,  and  other  optical  parameters  are  assumed  to  be 
uniform.  This  means  that  the  physical  thickness  of  a layer  would  be 
typically  too  large  to  permit  the  use  of  equation  (18)  for  repre- 
senting its  self-emission  levels.  The  apparent  dilemma  posed  by 
these  contrasting  requirements  can  be  resolved  by  using  specialized 
forms  of  equations  (22)  and  (23)  which  lead  to  the  definition  of  an 
emission  operator,  E,  analogous  to  S and  T.  The  existence  of  the 
operator  E,  which  is  not  internally  dependent  upon  the  temperature 
assigned  to  a layer,  allows  the  use  of  the  doubling  equations  to 
establish  S,  T,  and  E matrices  for  a layer  of  reasonable  thickness. 
Those  matrices  can  then  be  stored,  for  example,  on  magnetic  tape, 
and  never  need  to  be  recomputed  for  a given  cloud  model.  One  can 
therefore  define  a "basic  layer"  via  doubling,  and  use  its  S,  T,  and 
E operators  to  examine  the  effects  of  larger  scale  inhomogeneltles 
on  thermal  emissions. 

The  procedure  described  above  may,  of  course,  be  used  to  deal  with 
identical  layers,  in  which  case  the  subscripts  are  irrelevant.  If 
the  layers  are  homogeneous,  then  the  superscript  (*)  is  also  immaterial 

and  may  be  dropped.  Under  such  assumptions  I^t,  “ iti,,  and  one  may 

simply  write  I in  their  places,  yielding  [from  eqs.  (22)  and  (23)] 

ij  - [1  + TQ(1  + S)]  IE  s GIe  , (24) 

where  G is  a new  operator  given  by  G = 1 + TQ(1  + S).  It  will  then 
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also  be  evident  that  I, 


I*  for  an  Isolated  emitting  layer. 


Equation  (24)  provides  a powerful  tool  for  dealing  with  self-emission 
through  stacking  by  eliminating  any  need  to  determine  the  actual  values 

It  will  be  seen  that 


+ 4- 

of  I or  1 during  the  early  stages  of  doubling 
EE 


G is  formed  from  the  S and  T matrices  of  the  two  layers  being  combined 
in  a certain  step.  In  the  next  doubling  step  a new  G,  say  G'  will  be 
formed  using  the  combined  S and  T operators  which  may  be  called  S'  and 
T' 


Let  I be  the  result  of  the  first  doubling  step,  and  I'  be  that 
E £ 


of  the  next  step.  One  has,  then. 


GI. 


E’ 


and 


* 


G ' I_ 


(25a, b) 


Substitution  of  equation  (25a)  into  (25b)  then  yields 


El, 


where  E = G'G 


(26) 


The  new  operator  E thus  represents  an  "overall  emission  operator"  which 
operates  on  the  intensity  emitted  by  the  elementary  layer  used  at  the 
beginning.  As  successive  doublings  are  performed,  E may  be  updated 
through  multiplication  on  the  left  by  a G which  is,  itself,  fixed  by 
the  cumulative  S and  T operators. 


Like  S and  T,  E really  carries  several  parameters  and  is  a short  form 

for  its  true  representation  Ex (t ,y ') . Similarly  I « EI_,  actually 

A E E 


stands  for: 


-2-j  vt>vv 


(27) 


where  a.  is  the  Gauss-Legendre  weighting  coefficient  *or  abscissa  value 


V 


TREATMENT  OF  SURFACE  EMISSIONS  AND 
THE  CLEAR-AIR  LAYER  BELOW  A CLOUD 


This  section  will  show  the  ease  with  which  the  layer  stacking  method 
can  be  adapted  to  a variety  of  situations.  The  planetary  surface  is 
different  from  an  aerosol  layer  only  in  that  it  is  described  by  dif- 
ferent S,  T,  and  E matrices,  and  the  same  may  be  said  of  the  atmo- 
spheric layer  lying  between  the  cloud  and  the  surface  (or  above  the 
cloud) . 


The  surface  is  assumed  to  be  an  isotropic  emitter  and  scatterer  with 
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.ku etc  tcs+m* 


a reflectivity  and  an  emissivity  e ■ 1 - m . Such  a surface  can 
be  represented  by  a scattering  matrix  « WgU^/CAn),  and  a trans- 
mission matrix  Tjj  - 0. 

To  reasonable  accuracy,  a "clear"  atmospheric  layer  can  be  regarded 

as  neither  emitting  nor  scattering  appreciably,  but  as  simply  causing 

attenuation.  The  appropriate  S matrix  is  zero,  and  the  transmission 

matrix  T.  can  be  written  as 
A 


(t  ) 

VVij 


-T./U 


A'Mi 


ij* 


where  is  the  optical  depth  at  normal  incidence  of  the  atmospheric 
layer  at  the  wavelength  of  interest. 

By  use  of  the  layer  stacking  algorithm,  the  surface  and  clear-air 
layers  can  be  combined  into  a single  equivalent  layer  whose  trans- 
mission matrix  T^g  and  scattering  matrix  S ^ are  given  by 


T 


AS 


(28) 


and 


AS 


T ST 
A A 


(29) 


Let  l£  represent  the  emitted  intensity  at  the  planetary  surface  and 

i 

let  I*  and  X*  represent  the  upward  and  downward  intensities  at  the 
B B 

cloudbase.  With  these  definitions  one  can  then  write: 


li  = T 1+  + S 
B AS  S AS  B 


(30) 


and 


3 " 4 + Vb  * 


(31) 


where  Sc  is  the  reflection  operator  for  the  cloud.  Application  of  the 
Barkstrom-Querfeld  [A]  matrix  analysis  then  yields 


’ (1  - W1  + SAS^'  • 


(32a) 


and 


*£  “ ^ + SC(1  “ SAsV  SAS]IE  + SC(1  “ SAsV 


(32b) 
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Knowing  and  the  transmission  operator  T of  the  cloud,  it  is  a simple 

matter  to  see  that  the  portion  of  the  total  upward  intensity  emitted 
from  the  surface  which  is  seen  at  the  top  of  the  cloud  is  just 

[TC(1  ~ SASScrlTA]IS*  and  thSt  [TC(1  " SASScrlSAS]IE  tellS  h°W  “Uch 


of  the  upward  intensity  results  from  downward  emission  by  the  cloud, 
attenuation  by  the  atmospheric  layer  below  the  cloud,  and  subsequent 
upward  transmission  through  the  cloud.  In  fact,  the  latter  terms  also 
Include  all  multiple  reflections  between  the  cloud  and  the  surface. 

When  all  such  things  are  combined,  the  total  upward  intensity  seen 
when  looking  down  on  the  cloud  from  above  is  given  by 


up 


v1  - slssc>'Vte + [■ 


V1  - SASSC)'ll*]IS-  <33*> 


I = 1 + 
down  B* 


(33b) 


wherein  I*  and  1^  are  the  upward  and  downward  self-emission  intensities 
£ £ 

of  the  existing  cloud  in  the  absence  of  the  clear  air  and  the  surface 
Itself. 


VERIFICATION  OF  COMPUTER  CODE  CLEM 

The  previous  formulations  have  been  incorporated  into  a computer  code 
that  is  capable  of  predicting  theoretical  intensities  arising  from  cloud 
emissions  and/or  surface  emissions  passing  through  an  atmosphere  lying 
between  the  surface  and  cloud.  CLEM  is  described  fully  in  the  appendix. 

Verification  of  the  general  soundness  of  the  thermal  emission  code  has 
been  attempted  by  comparing  the  results  of  the  new  codes  with  some  of 
the  published  results  of  Yamamoto,  Tanaka  and  Kamltani  [5],  hereafter 
referred  to  as  YTK,  at  X ■ 10pm.  Since  the  particle  distribution  used 
by  YTK  appeared  to  be  the  same  as  Deirmendjlan' s [6]  cumulus  model 
C.l,  that  model  was  used  to  calculate  the  phase  function.  Table  1 
presents  the  results  of  one  such  comparison,  but  it  must  be  realized 
that  the  numerical  values  cited  as  representing  the  results  of  YTK  were 
taken  from  small  published  graphs  and  may  not  be  closer  than  5 percent 
to  their  actual  results.  The  agreement  appears  to  be  quite  satisfactory 
and  suggests  that  the  doubling  codes  contain  no  major  errors. 

Results  of  assuming  that  thick  clouds  can  be  regarded  as  isothermal 
(the  works  of  Zdunkowskl  and  Choronenko  [7]  and  of  YTK  [5],  for  example) 
have  been  examined  through  a series  of  computations  made  at  X ■ 11.5pm 
and  using  the  C.l  cloud  model  phase  function,  a variety  of  layer  thick- 
nesses, and  a lapse  rate  of  5*C  per  kilometer.  Each  layer  was  assigned 
the  temperature  that  would  prevail  at  the  center  of  the  layer  (the 
average  temperature).  In  addition  several  calculations  also  were  made 
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TABLE  1.  COMPARISON  OF  RESULTS  FROM  CLEM  WITH  THOSE  OF 
YTK  WITH  A SURFACE  TEMPERATURE  OF  30«c  AND  A 
CLOUD  TEMPERATURE  OF  -30°C  (A  - 10pm) 


Optical  Depth 

0.1 

1.0 

10 


zenith 
angle (°) 

0 

70 

0 

70 

0 

70 


Self-Emission 
Intensities 
(W  m~2  str-1  cm) 
Present 
Calculations 


1.2 

3.8 

11.1 

21.6 

31 

31 


YTK 

1.3 

5 

12 

13 

29 

29 


Total 

Upward  Intensities 

(W  m“2  str-1  cm) 
Present 
Calculations 


101 

89 

79 

A8 

31 

31 


YTK 

101 

91 

80 

46 

29 

29 


using  the  isothermal  assumption  and  taking  the  cloud  temperature  to  be  ) 

that  of  the  top  of  the  cloud.  Table  2 presents  the  results  obtained 
uftder  the  assumptions  that  surface  albedo  **  0,  that  cloud  base  tem- 
perature is  8°C,  that  surface  temperature  is  30°C,  and  that  absorption 
by  the  atmosphere  below  the  cloud  base  is  negligible.  In  the  table,  the 
quantities  listed  as  the  "apparent  temperatures"  are  the  temperatures 
which  would  be  inferred  if  the  upwardly  emitted  intensities  were  assumed 
to  be  produced  by  a blackbody.  The  data  shown  in  table  2 provide  a 
guide  as  to  the  fineness  needed  in  subdividing  the  layers  of  a cloud 
model,  and  also  indicate  that  the  assumption  of  isothermality  (with  the 
cloud  top  temperature)  does  not  lead  to  much  distortion  of  the  results. 

As  expected,  if  the  clouds  are  fairly  thin,  a substantial  portion  of  the 
upward  radiation  arises  from  the  surface,  making  the  apparent  tempera- 
ture quite  a bit  warmer  than  the  cloud  temperature.  Therefore,  both  the 
"graded  temperature"  cloud  models  and  the  isothermal  models  give  about 
the  same  value  (the  total  transmission  matrix  will  be  the  same  whether 
the  cloud  is  isothermal  or  not).  As  the  cloud  becomes  thicker,  the 
discrepancy  between  the  two  temperature  models  increases,  but  it  reaches 
only  about  1°C  for  a cloud  thickness  of  1200  m.  Also,  results  of  using 
step  sizes  of  100  to  200  m yield  results  very  close  to  those  obtained 
using  a 25-m  layer  thickness.  This  observation  suggests  that  even  if 
one  wants  to  use  a noniso thermal  model,  a layer  thickness  of  100  or  200 
m is  probably  small  enough. 

In  general,  as  the  clouds  reach  optical  depths  of  about  10,  the  isother- 
mal model  will  Vive  an  apparent  temperature  which  is  slightly  colder 
than  the  cloud  top  temperature,  while  a more  realistic  temperature 
profile  yields  a value  somewhat  (~0.6°C)  warmer  than  the  summit  tempera- 
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TABLE  2.  COMPARISON  OF  APPARENT  TEMPERATURES  FOR  ISOTHERMAL 
AND  NONISOTHERMAL  CLOUD  MODELS. 

Apparent  Temperatures* 

(Lapse  rate  - 5°C  per  km) 

Isothermal 

Computation  Step  Thickness  (m)  Single  Step  to 


Cloud 

Depth 

(m) 

Cloud-Top 

Temp. 

CO 

Optical 

Depth 

25 

50 

100 

200 

300 

the  Given  Depth 

T - T 

top 

25 

7.875 

0.249 

26.7 

— 

— 

— 

— 

50 

7.75 

0.498 

23.8 

23.8 

— 

— 

— 

— 

75 

7.625 

0.748 

21.3 

— 

— 

— 

— 

— 

100 

7.5 

0.997 

19.1 

19.1 

19.1 

— 

— 

19.0 

125 

7.375 

1.246 

17.3 

— 

— 

— 

— 

— 

150 

7.25 

1.495 

15.6 

15.6 

— 

— 

— 

— 

175 

7.125 

1.744 

14.2 





___ 

__ 

200 

7.0 

1.993 

13.0 

13.0 

13.0 

13.1 

— 

12.7 

300 

6.50 

2.990 

9.7 

9.7 

9.7 

— 

9.9 

— 

400 

6.00 

3.988 

7.8 

7.8 

7.8 

7.9 

— 

7.2 

500 

5.50 

4.984 

6.6 

6.6 

6.6 

— 

— 

— 

600 

5.0 

5.980 

— 

5.8 

5.8 

5.9 

6.0 

— 

800 

4.0 

7.974 

— 

4.6 

4.7 

4.7 

— 

3.9 

900 

3.50 

8.970 

— 

4.1 

4.1 

— 

4.3 

— 

1,200 

2.0 

11.961 

— 

— 

2.6 

2.7 

2.8 

1.8 

1,500 

0.5 

14.951 

— 

— 

1.1 

— 

1.3 

— 

1,800 

-1.0 

17.941 

— 

— 

-0.4 

-0.4 

-0.2 

— 

♦Surface 

temperature  +30*C,  cloud 

base  temperature  ■ +8 

°C;  surface  albedo  - 

o, 

X ■ 11.5pm 
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ture.  In  either  case,  however,  the  differences  between  the  actual  cloud 
top  temperatures  and  the  apparent  temperatures  are  again  so  small  as  to 
be  nearly  insignificant  if  t ^ 10. 

One  of  the  major  problems  to  be  encountered  in  temperature  retrieval  in 
future  automated  analyses  is  that  of  determining  whether  any  clouds  are 
present  and,  if  so,  finding  a way  to  estimate  their  optical  depth. 

These  pieces  of  information  are  necessary  in  assessing  the  relative 
contributions  to  the  observed  radiance  which  arises  within  or  below  a 
cloud.  It  is  conceivable  that  some  relevant  information  might  be 
obtained  by  making  the  radiance  measurements  at  more  than  one  infrared 
wavelength.  To  see  how  that  might  be  done,  let  the  earth's  surface  be 
assumed  to  be  a blackbody.  If  no  intervening  clouds  (or  only  very  thin 
clouds)  exist,  the  inversion  of  the  Planck  function  for  I,  and  I. 

A1  A2 

should  give  the  same  apparent  temperature.  The  same  is  true  if  the 
clouds  are  very  thick  since  their  behavior  does  become  similar  to  that 
of  a blackbody.  If  clouds  of  "intermediate"  depth  are  present,  the 
situation  may  be  different  since  the  transmission  and  emissive  proper- 
ties of  the  cloud  are  wavelength  dependent.  This  line  of  thought  has 
been  examined  by  comparing  apparent  temperatures  at  wavelengths 
A = 10pm  and  A^  = 11.5pm  for  several  cloud  models,  thicknesses,  and 

surface  temperatures.  A few  of  the  results  are  summarized  in  figure  3. 

A comparison  of  apparent  temperatures  inferred  at  wavelengths  of  10pm 
and  11.5pm  can  provide  information  of  the  type  suggested  above  if  the 
cloud  base  and  surface  temperatures  differ  appreciably  and  if  suffi- 
ciently precise  intensity  measurements  can  be  made. 

SIMULATED  SMOKE  STUDIES 

Theoretical  studies  have  been  made  with  thermal  emission  code  CLEM  to 
determine  the  flux  received  by  an  idealized  detector  of  variable  aperture. 
The  studies  were  made  to  determine  the  effects  of  surface  radiation 
passing  through  a simulated  cloud  produced  by  white  phosphorous  munitions. 
The  surface  temperature  was  set  at  30°C  and  50°C  while  the  cloud  tempera- 
ture was  held  constant  at  100°C.*  The  wavelengths  chosen  were  1.06pm 
and  10.6pm,  and  the  flux  into  a detector  of  aperture  180  degrees  and 
8 degrees  was  computed.  The  power  into  the  detector  may  easily  be  com- 
puted by  multiplying  the  flux  by  the  appropriate  surface  area  of  the 
detector. 

The  results  should  be  used  with  care  for  the  following  reasons: 

1.  CLEM  does  not  include  a "direct  beam"  but  rather  is  concerned 
with  multiple  scattering  effects  arising  from  thermal  emissions;  there- 


*This  temperature  was  for  white  phosphorous  smoke,  perhaps,  unrealis- 
tically higher  than  that  of  the  environment,  but  that  observation 
does  not  Impair  conclusions  related  to  the  way  in  which  a smoke  of 
such  a high  temperature  would  behave. 
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Fig.  3.  Differences  in  apparent  temperatures  at  two 
distinct  wavelengths  as  a function  of  cloud 
depth.  The  solid  curve  represents  a surface 
temperature  of  30#C  and  a cloud  base  temper- 
ature of  0*C;  for  the  dashed  curve  T =30#C, 
T£B=8°C;  for  the  dotted  curve  Ts=20°C, 
Tcb*0*C.  The  surface  reflectivity  was  zero, 
ana  the  cloud  model  was  Deirmendjian's  C.l. 


fore,  the  detector  is  only  "seeing"  the  thermal  radiation  that  is  emitted 
from  the  surface  and/or  cloud. 

2.  The  particle  distribution  [8]  may  vary  with  time. 

3.  In  lieu  of  reliable  values  for  the  refractive  indices  for  white 
phosphorous  smoke,  the  refractive  indices  of  sulfuric  acid  were  used  [9],. 

4.  The  detector  in  CLEM  is  idealized;  i.e.,  reflections  off  the 
detector  walls  and  internal  noise  are  not  accounted  for. 

Figures  4 and  5 (a  and  b)  present  the  ratio  of  surface  emission  to  sur- 
face emission  plus  cloud  emission  at  wavelengths  of  1.06pm  and  10.6pm 

3 _i 

for  a reference  density  of  1.2763  x 10  particles  cm  and  at  a surface 
temperature  of  30°C  and  50°C.  The  aforementioned  density  is  based  on 
values  from  Gomez  [8].  The  bottom  abscissas  (thickness  in  kilometers) 
of  figures  4a  and  b (or  figures  5a  and  b)  as  may  be  scaled  for  other 
densities  by 


N 


L (km)  = L (km)  -M.  = 1.2763  x 10 
New  Ref  N 

Npu 


-.3  LRef 


N 


New 


(34) 


Where  N is  the  new  density  in  par ticles/cm^ , and  L„  is  the  thick- 
New  Ref 

ness  in  km,  read  from  figure  4 or  5. 

Thus,  if  it  is  assumed  that  an  M15  (STD.C)  grenade  (0.9  lb  fill  weight) 

U 

is  used,  corresponding  to  5.203  x 10  particles,  and  that  the  initial 

3 -3 

burst  covers  8 m , the  average  density  will  be  65.04  particles  cm 

To  achieve  the  same  attenuation  as  the  reference  density  above  for 

t = 0.39(A  = 1.06pm,  T = 30°),  the  thickness  must  be  L,  (km)~9.8. 

s New 

Figures  6 and  7 present  the  flux  density  due  to  surface  emission  that 
has  penetrated  through  the  cloud  into  a 180-degree  and  8-degree  aperture 
detector  versus  thickness.  The  flux  density  due  to  cloud  emission  may 
be  found  by  using 


CF  ■=  SF/R  - SF,  (35) 

where  CF  is  the  cloud  flux  density,  SF  is  the  surface  flux  density 
(from  figure  6 or  7) , and  R is  the  ratio  surface/ (surface  + cloud) 
which  may  be  found  from  either  figure  4 or  5.  One  should  take  care 
that  the  same  detector  aperture,  surface  temperature,  wavelength,  and 
physical  depth  are  used  when  computing  the  cloud  flux  density  by  the 
above  formula.  Figures  6 and  7 may  also  be  scaled  in  the  same  fashion 
as  previously  mentioned  for  the  ratio  of  surface/(surface  + cloud) 
emission. 


Optical  Depth  (1.06um/10.6pm) 


Optical  Depth  (1 .06pm/ 10. 6pm) 
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Optical  Depth  (I .06um/10.6gm) 


ower  curve  corresponds  to  a detector  aperture  of  180  . The  surface  temp- 
30°C;  the  cloud  temperature  was  100°C. 


r-  ro 
oo  r- 


Thickness  (km) 

Fig.  5b.  Ratio  of  surface/ (surface+cloud)  emission  for  wavelengths  of  1.06ym  (dashed  lines)  and 
10.6ym  (solid  lines).  The  upper  curve  of  each  set  corresponds  to  a detector  aperture 
of  8#;  the  lower  curve  corresponds  to  a detector  aperture  of  180°.  The  surface  tenp- 
erature  was  30*C;  the  cloud  temperature  was  100°C. 


Optical  Depth 
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X=1 .06vim 


Thickness  (km) 

Log  surface  flux  versus  thickness  for  a wavelength  of  1.06ym. 
The  solid  curves  correspond  to  a surface  temperature  of  50#C; 
the  dashed  curves  correspond  to  a surface  temperature  of  30*C 
The  upper  curves  in  each  set  are  for  a detector  aperture  of 
18C^;  the  lower  curves  are  for  an  aperture  of  8*.  The 
cloud  temperature  was  100*C. 


Optical  Depths 
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a=  10.6pm 


Thickness  (km) 

Log  surface  flux  versus  thickness  for  a wavelength  of  10.6pm. 
The  solid  curves  correspond  to  a surface  temperature  of  SO® 
the  dashed  curves  correspond  to  a surface  temperature  of  30 
The  upper  curves  in  each  set  are  for  a detector  aperture  of 
180*;  the  lower  curves  are  for  an  aperture  of  8®.  The  cloud 
temperature  was  100®C. 
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The  first  point  that  should  be  noted  is  that,  for  a given  cloud  tem- 
perature, as  the  surface  temperature  decreases  the  probability  of 
detecting  the  surface  at  intermediate  optical  depths  decreases,  since 
as  the  temperature  of  a body  drops  so  does  its  resultant  thermal 
emission.  Secondly,  at  very  small  optical  depths  the  surface  emission 
is  recovered,  as  the  ability  of  the  cloud  to  absorb  or  emit  radiation 
is  severely  limited  by  its  "thinness."  Conversely,  at  large  optical 
depths  the  cloud  emission  is  recovered  as  the  surface  emission  is 
absorbed  by  the  cloud.  Clearly  these  results  are-  wavelength  dependent. 
Also  note  that  as  the  cloud  temperature  Increases,  ability  to  detect 
the  surface  emission  at  any  wavelength  will  be  decreased. 

In  examining  the  effects  of  the  detector  aperture,  two  limiting  cases 
are  considered:  (1)  small  optical  depth,  and  (2)  infinite  optical  depth. 
Regardless  of  the  detector  aperture,  at  small  optical  depths  absorption 
and  emission  of  radiation  by  the  cloud  will  be  practically  negligible 
due  to  the  "thinness"  of  the  cloud.  Thus,  the  surface  radiation  will 
pass  through  this  thin  cloud  with  practically  no  attenuation,  resulting 
in  a surface/ (surface  + cloud)  ratio  that  is  nearly  unity.  For  an 
optically  thick  medium,  the  cloud  emission  is  an  overriding  factor  at 
all  angles  regardless  of  aperture  size;  as  the  cloud  becomes  optically 
thick  not  only  does  the  cloud  emission  approach  a uniform  value  (black- 
body  radiation)  but  it  will  also  absorb  any  radiation  incident  from 
below.  However,  the  first  of  these  effects  is  dominant. 

At  intermediate  optical  depths  the  situation  is  more  complicated:  the 
cloud  emission  will  be  less  at  smaller  angles  due  to  a shorter  path 
length  and  the  surface  emission  seen  will  be  greater  at  small  angles 
due  to  less  material  lying  along  the  path  length  (less  absorption  and 
scattering).  Thus  at  intermediate  optical  depths  a smaller  detector 
aperture  will  "see"  (1)  effectively  a larger  surface  component  and  (2) 
a lesser  cloud  component  than  a larger  detector  aperture.  This  result 
is  dependent  on  the  wavelengths  and  the  index  of  refraction,  n.  For 

example  at  1.06pm,  n = 1.426  - 1.4  x 10"^i  and  at  10.6pm,  n = 1.953  - 

0. 4681  for  figures  4 and  5.  Therefore,  the  absorption  at  1.06pm  is 
much  less  than  at  10.6pm  and  the  results  agree  with  the  above  analysis, 

1. e.,  at  intermediate  optical  depth  the  results  of  the  8-degree  detector 
agree  with  the  180-degree  detector  much  better  at  1.06pm  than  at  10.6pm, 
due  to  the  difference  in  absorption  at  these  respective  wavelengths. 
However,  one  should  keep  in  mind  that  although  the  ratios  of  surface/ 
(surface  + cloud)  emission  may  be  approximately  the  same  for  the  8- 
degree  detector  and  the  180-degree  detector,  for  any  optical  depth,  the 
actual  flux  may  differ  by  orders  of  magnitude,  as  may  be  seen  in  figures 
6 and  7. 

In  summary,  at  small  optical  depths,  the  surface  may  be  detected  with 
little  problem  at  any  wavelength,  while  at  large  optical  depth  there 
is  little  chance  of  detecting  emission  from  the  surface.  At  intermediate 
optical  depth,  the  ability  to  detect  surface  emission  is  highly  dependent 
upon  the  index  of  refraction  of  the  medium  under  consideration.  Also, 
the  thermal  emission  arising  from  the  cloud  must  be  accounted  for. 
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APPENDIX 


DESCRIPTION  OF  CLOUD  EMISSION  COMPUTER 
CODE  CLEM 

CLEM  is  a thermal  emission  routine  that  will  calculate  intensities, 

via  the  doubling  method,  at  the  top  and  bottom  of  a cloud  base  due  to 
surface  emission  and  cloud  emission.  Because  CLEM  was  constructed  to 
operate  at  Infrared  wavelengths,  no  direct  beam  is  included.  Treatment 
of  continuum  absorption  by  water  vapor  for  saturated  layers  is  avail- 
able in  the  program.  The  following  description  is  aimed  at  understand- 
ing how  the  program  operates,  what  arrays  are  used  for  what  variables, 
and  a general  description  of  the  computer  code.  It  should  be  noted 
that  layers  are  added  to  the  bottom  of  the  stack  instead  of  the  top 
within  the  program. 

Subroutine  EMGUST  is  called  to  determine  the  Gauss-Legendre  weights 
[W(I) ] and  abscissas  [EMU (I)]  by  Newton-Ralphson  iteration  for  the 
interval  (-1,1).  After  the  return  to  MAIN  the  abscissas  are  shifted 
into  the  Interval  0,1  and  SIG(I)  « 1/p^,  and  W(I)  = W(I)/ui,  where 

■ cos  9^  are  calculated  and  stored  along  with  the  angles  themselves 

[EMO(I)].  As  the  index  I increases,  9 will  go  from  approximately  90 
degrees  to  approximately  0 degrees  (u  = 0 1).  The  next  read  is: 

NLYRS,  the  number  of  layers  of  thickness  THICK  (in  kilometers);  ALBS, 
the  surface  albedo;  TMSC,  the  temperature  of  the  surface  in  degrees 
Centigrade;  MODEL:  0 -*■  skip  subroutine  VAPOR,  1 ■+  AFCRL  tropical  atmo- 
spheric model,  2 AFCRL  midlatitude  summer  model,  3 AFCRL  midlati- 
tude winter  model,  4 -*-1962  US  standard  atmosphere  model,  99  cloud 
saturation  only;  HEIGHT,  the  cloud  base  height  in  kilometers;  PUNCH, 
if  less  than  one,  punch  output  for  the  plotter;  TAUR(I)  for  I ■ 1, 

NLYRS,  temperature  in  degrees  Centigrade  of  layer  number  I.  This  pre- 
vious data  must  be  user  supplied,  whereas  all  subsequent  data  are  pro- 
duced by  the  Mle  program.  The  next  read  is:  NWAVE,  the  number  of 
wavelengths  under  consideration;  NN,  the  number  of  Legendre  coefficients 
GNU,  the  wavenumber  under  consideration;  and  CT,  the  extinction  coef- 
ficient. 

The  next  read  is  for  the  Legendre  coefficients  [0L(I)],  where  I - 1,NN 
and  NN  equals  the  number  of  coefficients  (which  cannot  exceed  32) . 

The  phase  function  should  be  normalized  so  that  the  first  Legendre  co- 
efficient 0L(1)  is  equal  to  the  single-scattering  albedo  uTq.  The  next 

call  is  to  VAPOR,  which  will  set  up  the  proper  atmosphere  (below  the 
cloud  base)  according  to  whether  MODEL  is  equal  to  1,  2,  3,  or  4,  and 
return  TAUA  the  optical  depth  of  the  atmosphere  below  the  cloud  base 
adjusted  for  water  vapor  absorption  at  temperature.  Next  compute  TMPS, 
the  Kelvin  temperature  of  the  surface  and  TAUI,  the  optical  thickness 
of  the  basic  layer  (cloud).  Then  a few  manipulations  are  performed  to 
determine  th^  number  of  doublings  (NTI)  necessary  to  get  to  the  final 
optical  depth.* 
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Subroutine  CLSTK  is  a doubling  routine,  which  basically  will  double  up 
to  the  basic  cloud  (or  layer)  thickness  and  determine  the  total  (for 
this  layer)  scattering,  transmission,  and  the  multiple  scattering  opera- 
tors for  thermal  emissions.  Note  that  CLSTK  is  not  concerned  with  the 
atmosphere  or  surface,  but  only  with  the  cloud  (basic  layer).  Unless 
otherwise  noted,  here  and  elsewhere  the  incident  angles  are  labeled  J 
and  the  exit  angles  are  I.  Computations  are  performed  to  determine  the 
Legendre  polynomials  [PL(I,J)];  find  the  beginning  layer  thickness  [TAU]; 
find  the  exponential  attenuation  at  TAU  for  a given  incident  angle 
[ETS(I)];  set  F(I,J,1)  equal  to  the  identity  matrix;  find  the  phase 
function  needed  for  initial  computation  of  the  scattering  matrix  [nega- 
tive argument,  ends  up  in  P(I,J)]  and  for  the  transmission  matrix 
[positive  argument,  ends  up  in  R(I,J)];  finally  calculate  the  initial 
scattering  and  transmission  matrices  [S(I,J,1)  and  T(I,J,1)]  and  add 
in  the  exponential  term  on  the  diagonal  to  T(I,J,1).  As  NEWOLD  = <p , 
in  this  call,  now  proceed  to  double  the  appropriate  matrices. 

The  basic  doubling  loop  for  CLSTK  is  loop  26,  and  the  general  format  is 
to  calculate  the  appropriate  matrices  from  the  previous  ones.  This  cal- 
culation is  accomplished  by  indexing  S,  T,  and  F (the  scattering,  trans- 
mission and  "thermal"  operators)  with  a third  subscript,  equal  to  1 or 
2.  On  the  first  pass  LA  = 1,  LB  = 2;  and  on  subsequent  passes  their 
order  is  interchanged,  as  the  "future"  operator  is  calculated  from  the 
previous  one. 


Thus  the  first  operation  that  is  performed  is:  Q * 1 - S^S^.  The  Q 

matrix  is  then  inverted  by  either  (1)  taking  the  first  two  terms  of  a 
MacLaurln  expansion  (if  Q is  small  and  IFAST  has  been  set  # 0) , or 
(2)  calling  a Gauss-Jordan  matrix  inversion  routine  [EMINV].  In  either 

case  Q is  returned  as  (1  - S.^S2)“^.  Starting  now  at  label  14  compute 

A(I,J)  = (1  - s^S2^_1ti  anci  = TiS2*  The  aPProPriate  matrices 

are  now  available  to  determine  the  "new,"  or  doubled,  values  for  S and 
T.  Thus  (note  the  third  subscript  here)  T(I,J,2)  * T2(l  - 

the  "new"  value  of  T,  and  S(I,J,2)  = + T^S^l  - S.^)”1^,  t*ie  "new" 

value  of  S.  The  rest  of  this  do  loop  (26)  is  used  to  compute  the  "new" 
thermal  emission  operator  F.  Thus,  first  compute  C(I,J)  = 1 + S,  where 
S is  the  original  scattering  matrix,  not  the  "new"  or  doubled  value. 

At  this  point,  also  drop  the  subscripts  on  T and  S (i.e.,  T^  » T, 

■ S,  etc.);  as  for  a homogeneous  layer,  thermal  emission  is  iso- 
tropic. Next  compute  Q1(I,J)  - (1  - SS)"^(1  + S)  and  then  G(I,J)  * 

1 + T(1  - SS)_1(1  + S).  The  final  part  of  this  loop  (26)  computes  the 
new  or  doubled  value  of  the  thermal  emission  operator:  F(1,J^2)  - 

[1  + T(1  - SS)_1(1  + S)]  x F(I,J,1),  where  one  should  recall  that  on 


29 


I 


r 


T 


the  first  pass  F(I,J,1)  = 1.  The  optical  depth  is  then  incremented 
and  the  doubling  concinues  until  the  final  optical  depth  of  the  basic 
layer  is  reached. 

After  the  doubling  loop  is  complete,  the  final  thermal  emission  matrix 
operator  for  a homogeneous  layer  is  now  stored  in  F;  S and  T still  have 
their  usual  cleanings.  Subsequently  the  subroutine  will  return  values 
to  MAIN  for  the  final  (doubled)  basic  layer  via  a switch  and  equivalence 
statement.  Thus  the  returned  values  for  the  basic  layer  are:  A -*■ 
current  "doubled"  value  of  S;  B -*■  current  "doubled"  value  of  T;  and 
P -*■  current  "doubled"  value  of  F.  The  switch  determines  that  the  cor- 
rect values,  and  not  the  previous  values,  are  placed  in  the  proper 
arrays. 

After  the  return  to  MAIN  the  albedo  for  single  scattering  for  the  cloud 
is  saved  in  ALBC.  Now  if  MODEL  = 0 VAPOR  is  effectively  bypassed.  If 
not,  set  TMP  equal  to  the  Kelvin  temperature  of  the  cloud  and  call  VAPOR, 
which  (last  argument  = 99)  will  change  the  physical  depth  of  the  cloud 
(rather  than  the  optical  depth)  to  account  for  absorption  by  water  vapor 
molecules  at  100  percent  relative  humidity.  TMP  is  returned  as  an  addi- 
tional extinction  coefficient  which  is  used  to  reduce  the  thickness 
[TAUT]  of  the  basic  layer. 


Now  0L(I)  and  Q1(I,1)  are  computed  as  the  thermal  emission  due  to  a black- 
body  whose  emission  is  1 and  albedo  of  single  scattering  is  0.  (Recall 
that  here  P = F.)  The  surface  emission  (ESURF)  is  then  calculated  by 
calling  PLANCK,  which  for  a given  wavenumber  (GNU)  and  temperature  (TMPS) 
will  return  the  blackbody  emission  as  BNUS.  Note  that  at  the  ground  sur- 
face tau  is  considered  infinite.  The  blackbody  radiation  of  a cloud  of 
infinite  optical  thickness  (ECLOUD)  is  then  found  by  the  same  procedure 
as  described  for  the  surface  emission,  but  at  the  temperature  for  our 
first  basic  slab  [TAUR(l)  + 273.16].  The  next  loop  (28)  first  finds  the 
interpolated  thermal  intensities  of  an  infinitely  thick  cloud  whose  black- 
body emission  is  1,  for  0 degrees  -*■  75  degrees  in  5-degree  increments. 

This  quantity  is  returned  from  the  Interpolation  routine  [EMNTRP]  as 
ALBUP(l).  The  true  thermal  emission  due  to  the  cloud  at  the  appropri- 
ate temperature  and  optical  depth  is  calculated  and  stored  in  Q1(I,5), 

ETS(I),  and  Q1(I,6).  Note  that  here  OL  - F* (1  - e-t/y)  and  ECLOUD  - 
(1  - d>o  )Bv(Tc).  Within  the  same  loop  the  combined  scattering  matrix 
c 

for  transmission  (downward)  by  the  atmosphere,  scattering  by  the  surface, 

and  retransmission  (upward)  by  the  atmosphere  as  seen  at  cloud  base 

(S  ) is  calculated  and  stored  in  PL(I,J).  Subsequent  to  this,  but  still 
AS 

within  the  same  loop,  there  is  an  array  transfer;  matrices  G,  BS,  and 
AS  are  set  equal  to  A (combined  value  of  S,  the  scattering  matrix  for 
the  basic  layer);  matrices  AT,  BT,  and  Q2  are  set  equal  to  B (combined 
value  of  T,  the  transmission  operator  for  the  basic  layer). 


Next,  subroutine  EMPR60  is  called,  which  essentially  processes  data, 


i 

] 


J 

I 

j 
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includes  surface  effects,  and  writes  output  for  the  first  layer,  in- 
cluding surface  and  atmospheric  effects.  First,  intensity  due  to  the 

surface  as  seen  at  the  cloud  base  is  calculated  tEYUP(I)  = (1  - w ) 

°S 

B(Tg)T^] , and  the  matrix  Q becomes  (1  - SASSC)»  which  is  then  inverted 

by  EMINV  and  Q is  returned  as  (1  - S S )-^.  New  matrices  B = 

AS  C 

T (1  - S S ) ^ and  P = S (1  - S S ) ^ are  then  calculated,  and  sub- 
L Ab  C C Ab  C 

sequently  the  matrix  R - T (1  - sACsr)_1S  is  also  computed.  At  this 

L Ao  L Ab 

point,  note  that  if  u “0  then  the  scattering  operator  goes  to  zero 

s 

and  thus  so  does  S , the  combined  transmission-scattering-retransmission 
Ab 

operator  for  the  atmosphere  and  surface.  In  this  case  P B + T^, 

R ->  0,  and  FL  -*  0,  where  an  appropriate  switch  has  been  incorporated  in 
the  program.  The  subscript  c represents  the  cloud.  Next  Q1(I,7)  = 

T (1  - S S ) T.l|,  the  upward  Intensity  due  to  the  surface,  looking 

C Ab  C A b 

down  on  the  cloud  from  above  is  found.  Note  that  EYUP  * T I*  ®lready  has 

A b 

the  atmospheric  attenuation  included.  Q1(I,11)  ■ T (1  - S S )-^S4  il, 

^ C AS  c AS  E 

+ I , the  downward  thermal  emission  due  to  the  cloud,  which  has  been 
b 

multiply  reflected  and  transmitted  by  the  cloud,  atmosphere.,  and  surface 

is  computed.  Then  the  matrix  R ■ S (1  - S S )_1S  is  calculated. 

C Ab  L Ab 

The  final  intensities  at  the  Gauss-Legendre  angles  are  then  calculated: 

Q1(I,10)  - PxEYUP  - s (i-S  £5  )_1T  I*  : the  surface  emission  in  the 

C AS  C AS  downward  direction  as  seen 

at  the  cloud  base 

Q1(1, 8)  - RxQ1c+Ql1;*s(n(l-sACSr)-lsAcIp+Iv  : the  cloud  emission  in  the 
5 5 C A downward  direction  as  seen 

at  the  cloud  base 

Q1(I,9)  - Ql  «l  -(1-S  S J^S  T IXS  Ijl+lj  : the  total  emission 
8 10  AS  C C A S C AS  E E ln  the  downward 

direction  as  seen 
at  the  cloud  base 

Qiu,i2> . ^\^ca-s^c)-hAi^ca-s^c)-h^i^  , th. 

total  emission  ln 
the  upward  direc- 
tion as  seen  at 
the  cloud  top 
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and  also 


Q1  (1,13)  - Qly/ESURF  - Tc(1-SasSc)-1TaI+/[(1-(L08)Bv(Ts)] 

Q1(I, 14)  - T 

Apparent. 

In  summary  then  the  Q1  matrices  are  as  follows: 


Ql^  -*■  upward  Intensity  at  TOC  due  to  the  surface  emission. 

Ql^  ■+  upward  Intensity  at  TOC  due  to  downward  emission  by  the  cloud. 

Ql^g  ■+  downward  Intensity  at  BOC  due  to  the  surface  emission. 

Qlg  -*■  the  downward  Intensity  at  BOC  due  to  multiply  scattered  cloud 
emission  plus  the  downward  cloud  emission. 


Qlg  -+  the  total  downward  Intensity  of  BOC. 


Q1^2  + the  total  upward  Intensity  at  TOC. 

The  last  part  of  this  routine  simply  Interpolates  [via  EMNTRP]  the  Ql's 
and  returns  them  as  ALBUP ' s , and  then  prints  the  necessary  quantities. 
There  Is  a one-to-one  correspondence  between  Q1  and  ALBUP  for  the  second 
Index  of  Ql. 


Now  return  to  MAIN  where  the  "stacking"  loop  is  entered  and  some  matrices 
are  rearranged:  Q1(I,3)  ■ Q1(I,4)  ■ thermal  emission  matrix  for  the 

next  layer  - G(1  - e"T/U)(l  -i  )By(Tc);  Q1(I,1)  - Q1(I,2),  the  thermal 

c 

emission  matrix  for  the  previous  layer;  A»P«S;B«R-T;AS-BS“G 
AT  • BT  ■ T.  Note  that  only  on  this  first  pass  does  A ■ P,  etc.  In  sub- 
sequent passes  j4  S2  and  therefore  A + P,  etc.  Next  find  the  total 

thickness  (kilometers)  of  the  "final"  layer  [TAUU],  and  after  calling 
VAPOR,  the  thickness  of  the  new  upper  layer  (kilometer)  [TAUL].  From 
here  go  to  subroutine  ADDER. 


ADDER  basically  will  combine  a "basic"  layer,  with  the  previous  layer 
(note  that  the  previous  layer  means  from  the  cloud  base  upward,  l.e., 
the  atmosphere  and  surface  are  not  Included  here,  but  rather  In  EMPR60) . 
ADDER  will  alao  combine  homogeneous  and  inhomogeneous  layers. 

The  first  loop  In  ADDER  essentially  rearranges  variables: 

EYUPU  - 1*^;  EYDNU  ■ 1*^;  EYUPL  ■ I*  ; EYDNL  - 1*^;  where  the  subscript 
1 denotes  the  previous  layer  and  2 denotes  the  upper,  or  new,  layer. 
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•k  * * * 

Also,  transfer  OS^;  D“T^;  E-S^;  F-T^;  G-S^;  H'T^:  note  that  when  ADDER 
is  called  A ■ and  B * T^.  Next  find  Q and  invert  it  so  that  Q ■ 

(1  - S*S2^~^*  Then  set  Q1  * Q and  R ■=  (1  - S*S'2)”^"T^.  Subsequently, 

P - S2(l  - S*S2)_1T1;  AT  - T2(l  - S*S2)"1T1  (which  is  the  "new"  T) ; and 

AS  - T*S2(1  - s*s2)~1T1  + S (which  is  the  "new"  S) . Now  if- INHOMO  - 2, 
then  the  layers  are  inhomogeneous  and  S*  ^ S^,  etc;  otherwise,  BS  ■ AS 
and  BT  ■ AT.  If  the  layers  are  inhomogeneous,  the  following  occurs: 

Q = (1  - S2S*)-1;  R = (1  - S2S*)'1T2;  BT  = T*(l  - S2S*)-1T2  (the  new 

T*);  P = S* (1  - S2S*)-1T*;  BS  - S2  + T2S*(1  - <the  new  s*) • 

From  here  the  rest  of  the  routine  is  concerned  with  finding  the  cloud 

intensities  for  the  combined  layers.  Thus,  A = S (1  - S*S  )-^;  B *= 

2 12 

T2(1  ~ S1S2)_1;  G = T1S2(1  " SiS2)_1;  E = T2(1  ' S1S2)_1S1*  N°W 
INHOMO  = 0 implies  that  the  two  layers  are  at  the  same  temperature  and 
thus  will  have  the  same  blackbody  emission.  If  this  is  not  so,.  H = 

T*S2(1  - S*S2)-1S*  + T*  is  calculated.  Next  the  resultant  thermal 

emission  intensity  in  the  downward  direction  is  found:  Q1(I,6)  - 

T2(l  - s*s2)  ls*I2E  + T2^  - SlS2^  1]C1E  + I2E*  and  lf  the  layers  are  at 
the  same  temperature  Q1(I,5)  = Q1(I,6);  else  Q1(I,5)  ■ T*S2(1  - S*S2)-1I^ 

+ [T  S (1  - S S )-1S  + T ]I,  +1.  , the  resultant  thermal  emission  in- 

1 Z 1 Z J.  1 Zb  Ih 

tensity  in  the  upward  direction.  Now  return  to  MAIN. 

The  return  to  MAIN  is  still  in  the  doubling  loop  and  processes  are  con- 
tinued until  the  final  optical  depth  is  reached.  The  program  will  loop 
back  to  502  CONTINUE  if  the  number  of  wavelengths  is  greater  than  one; 
in  this  case  new,  values  for  NN,  GNU,  CT,  and  the  Legendre  coefficients 
must  be  supplied  for  every  new  wavelength. 

INPUT  FOR  CLEM 

CARD  1 Input  Symbols 

NLYRS,  THICK,  ALBS,  TMSC 
[F0RMAT(I5,  3F10.6) ] 

NLYRS  - number  of  layers  of  thickness  THICK 
THICK  - thickness  of  individual  layers 
ALBS  - surface  albedo 
TMSC  - surface  temperature  (*C) 
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CARD  2 


Input  Symbols 


CARD  3 


CARD  4 


CARD  5 


CARD(S)  6 


MODEL,  HEIGHT,  PUNCH 
[FORMAT (15,  2F10.6) ] 

MODEL  ■ 0 for  no  atmospheric  continuum  absorption 

1 for  AFGL  Tropical  atmospheric  model 

2 for  AFGL  Midlatitude  summer  atmospheric 
model 

3 for  AFGL  Midlatitude  winter  atmospheric 
model 

4 for  US  1962  standard  atmospheric  model 
HEIGHT  - height  of  atmosphere  (£  9 km) 

PUNCH  ■ if  £ zero,  punch  "retrieved"  temperatures 
at  preselected  angles  after  each  doubling. 

Input  Symbols 

TAUR(I) , I - 1,  NLYRS 

[ FORMAT (10F8. 2)] 

TAUR  ( ) ■ true  temperature  of  individual  layers  (°C) 

Input  Symbols 


NWAVE 

[FORMAT (15)] 

NWAVE  * number  of  wavelengths  under  consideration 

Input  Symbols 

NN,  GNU,  CT 
[FORMAT(I5,2D12.6) ] 

NN  - number  of  Legendre  coefficients 

GNU  - wavenumber  (units  cm-*) 

CT  ■ extinction  coefficient  (km“^) 

Input  Symbols 


OL(I) 

[FORMAT  D25.14] 

OL(I)  is  the  Itl*  Legendre  expansion  coefficient  for 
the  phase  function. 

NOTE:  There  should  be  NN  cards  of  this  type. 

Subsequent  Cards:  NWAVE  additional  sets  of  type  5 and  type  6 cards  are 
required  if  NWAVE  > 1. 
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PROGRAM  LISTING 

CLEM 


1 

2 

3 

<1 

5 

6 
7 
n 

9 

in 
1 1 
12 
• 3 
II 
IS 

14 
17 

15 
l» 
20 
21 
22 
23 

29 

25 

26 
27 
?S 
27 

30 

31 

32 

33 
39 
35 
3* 

37 

38 
37 
SO 
S| 
S7 
S3 

ss 

S5 


COMPILER  I 0 1 A G ■ 3 ) 

IMPLICIT  REAL»B(A-H,o-7) 

C PROGRAM  CLEM70  • • ADAPTFO  FROM  CLEM6S ( 26  MAR  1977  I 17  OCT  1977 

C THIS  VFRSlON  CONCENTRATES  ALL  TEMP  DATA  IN  A FEW  CAROS. 

C ALSO  INCLUDES  H20-VAP0R  ARSORP.  PFLOW  ANO  WITHIN  THE  CIOUD 

C WILL  NOP  LOOP  ON  WAVELFNGTH  - RCS 

C UNITS  FOP  ECLOIID  AND  ALL  ’INTENSITIES'  ARE  MILLI  WATTS/ 

C ISO.  MFTER-SR-  INVFRSE  CENTIMETER! 

c 

1 FORMAT!  I5.3FI0.6, I I ) 

2 FORMAT  I /'  CONTINUUM  ABSORPTION  RV  WATER-VAPOR  WITHIN  THF  CLOUD  TIE 

• LOS  K I VAPOR  I • ’.IP0I2.S,’  PFR  KM. ’/I 

3 FORMAT!/’  LAYER  NO.  • , I 7 , 5 X , • TH I C K NESS  ■ ’ , F9 . K , S X , • T FMp  • ’,F7.3 

• , 3X  , ' BNUC  » ' • |PD| | • 9 | 3 X . ’FCLOUO  » ’.013. 6/1 

S  FORMAT!  ’ WAVENUMBFR  - ’,F  17. S,’  CM  -I.  SURFACE  ALBEDO  • ’, 

• IPDT3.6,’.  CLOUD  ALBFDO  FOR  SINGLE  SCATTFRING  « • .019.7/) 

5 FORMAT!  • StlRFACF  TFMP  * • , F R . 3 . 5 X , • RNUS  ■ • , I PQ  I I • 9 . 5 X , ’ESURF  • 
•’,DT3.6) 

6 FORMAT! /8I7X ■ IHL  *SX .SHOMFGA ,SX ) / I 

7 FORMAT  I R1 IX  , 12,F 13.9) ,9X ) 

R FORMAT! /••••  PHASE  FUNCTION  EXPANSION  TRUNCATED  TO  3’  TERMS  •••■) 

9  FORMATHH  ,’FOR  A TOTAL  DEPTH  OF  ’ »F8.3,’  KM  •/) 

10  F ORM  A T I | H | , • WAVENUMRER  ■ ’.GII.S,’  EXTINCTION  COEFFICIENT  ■ ’. 

♦Gll.Sf*  KM-I’I 

11  FORMAT!  ///’  OTHER  PARAMETERS  ’/) 

12  FORMAT!/’  INTENSITIES  FOR  THE  BASIC  SLAP  IF  ECLOUO  ■ I •/> 

13  FORMAT!’  ANGLF  ■ ’.F5.2,’  OEG.’.SX,’  INTEMS.  • ’,F9.S) 

1 7 FORMAT!’  THE  PHASE  FUNCTION  EXPANSION  CONTAINED  ’.13.’  TERMS  •) 

I R EORMATI IQFA.2) 

I 9 EORMAT 1075. | R I 

20  FORMAT! |HI  ,/) 

21  FORMAT!/! 

100  FORMAT  I 15) 

101  FORMAT  115, 2D  17.6) 

C 

REAL  TSS,fST,T,03,ISNGL 

DIMENSION  S(I6,|6,7>,T<|6,I6,2),A(|6,I6),B(|6,I6>,QI(|6,!6) 
DIMENSION  ASHA, 16),  BS  l|6,|6>, ATH4,|6),RTI|6, 16), ALP  UP  116) 
DIMENSION  nLl37)>EMUI|6),S|G(|6),WI|6),ETSI|6),EM0l|6l,Y(l6) 
DIMENSION  PL (16.16), PI|6,I6).RI|6,I6)>0(I6, |6I, APTII6) 

DIMENSION  G( 16, 16 ) , ANG! 16) ,07! |6  • 16  ) ,EYUP< I 6 I , E Y ON < I A ) , T AUR ( SO ) 

FOtll  VALENCE  IS!  | ) ,AS(  | ) ) , (SI257)  ,RS(  I ) I 

EQUIVALENCE  IT! I I ,AT ( | I I , (TI2571  ,RT  I I I I 

COMMON  S • T , A ,B  ,0  ,QI  ,P  ,R  ,PL  ,0L  ,EHtl,FMO.W  ,ETS  ,S  IG 

COMMON  FM,  TAIIU.EMUO.TAUL  , T AU I , T AlIT  , AL  BUP  . ECLOUD  , mm  , NN  , L . L L 


1 


R6 

97 

98 

99 
50 


COMMON  NT  I ,NTF,KO 

THIS  STMT  AND  THE  NEXT  WILL  NOT  APPEAR  |N  FURTHER  PG«  LISTINGS 

CALL  10 

NCOUNT-O 

L ■ I 6 


I 

I 
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Cl  EM 


*1 

PI*3»I9|59765 35 8980*00 

52 

FP-r.Onvno/ ( 9.0D+P0  *Pt) 

53 

CALL  E«GUST 

59 

no  7?  ! « 1 .L 

«S 

r Mo ( I ) * ( i .ooo+nu-EMo ( i ) ) / 7. on *00 

56 

ft  1!  ) * S T G ( I) 

57 

5 1 G ( I ) * I .Oo+oo/emu ( i i 

58 

W(|)=VMt)*5lG(T> 

59 

77 

r m o ( i ) = |80.on>ni. *dacos(emu(  i » i/pi 

60 

REAnl  05,|)  Ml.  Y95, THICK, AL85.TMSC 

61 

RE  AO » 05, 1 IMOOEL. HEIGHT, PUNCH 

6? 

I PIINCH*0 

63 

IE  ( PUNCH. LE  • 0 • ) | PUNf  H»  | 

69 

READ((J5,|8)  { T A HR ( I > ,I  = | , NL  Y R 5 ) 

65 

READ  <O5,l*l0>  NWAVE 

66 

502 

CONT I NUE 

67 

READ  (5, 1 0 | ) NN,GNU,CT 

68 

WRITE  f 6 • 1 0 ) GNU  , C T 

69 

WRITE(6,»7)NN 

70 

1 E ( NN  .I  E • 37  ) GO  TO  73 

▼I 

WR I TE ( A , 8 ) 

72 

N N * 3 7 

73 

23 

no  29  I * * t NN 

9*» 

79 

PEAD!u5,|9,eMD«29|)0l.  II) 

75 

291 

WR  | TF ( 6 , 6 ) 

76 

WR!TE(6,7) C I ,0L(  I ) *1*1  *NN) 

77 

WR  I TE  ( 6 , 1 1 ) 

78 

C 

79 

CALL  VAPOR! ME  I GHT, GNU,  TAMA,  MODEL* 

80 

TMP5  = TM5C*773.|6D^r»C 

8| 

TAM! *TH1 CK*f T 

82 

1 5NGL*  1 DLOG ( T AU I * / ( DL 0 G 1 7 • n *00 ) ) ) 

93 

TST*?0.0* 1 5NGL 

89 

NTP«IEIX(T5T) 

85 

RNTE*NTF 

86 

TST«TST-RNTE 

87 

IE(TST.GE«O.IO)NTE*NTE*l 

88 

NT  I =NTP 

09 

CALL  CL5TKIO,D) 

90 

M T E * I 

9 | 

A L 8 C * 01  ( 1 ) 

97 

1 E < MODEL »NE .0 ) GO  TO  795 

93 

TMP-0.0 

99 

GO  TO  796 

99 

795 

TMP  -TAUR ( 1 ) *273. 1 6 

96 

CALL  VAPORIHF I GHT , GNU , TMP , 99 ) 

97 

796 

T AU  I »T  AIIT 

98 

TMPO * T MP 

99 

TAUT*  TAIII /(CT*TMP) 

|Op 

DO  76  I * 1 *L 
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CLEM 


I 


•• 


) 

1 

i 

I 


101 

102 

103 
|P3 
105 

104 
107 
IPS 
10? 
I 10 
I I I 
112 

113 

no 

115 

114 
117 
II" 
II? 
I»0 

171 

172 

173 
17? 
125 

174 

177 

178 
I »? 

130 

131 
I 32 
133 
I 3? 
I 35 
I 36 
137 
I 18 
I 3? 
I?0 
l?l 
I ?? 
I 33 
l?3 
I "5 

136 

137 

138 
13? 
150 


T A «0 • 0*00 

00  25  J«|  ,L 

25  TA»T**5  IC,  ( J)  »W<  J » •EM1H  Jl  *P  I ! , J|  • (I  .d  + OO-ETS  ( J)  >/»(  I I 
Ql I t • I 1-TA 
76  OL ( t I *TA 

*fMTE<A,3lGNll,AL8S,ALBr  i 

CALL  PLANCK  ( GNII  .IMPS  ,BNUS  I 
FSURF- ( | .0*00-AL8S  t •BNIIS 
»R  t TF ( 6 ,5  I TMSC ,BNUS ,FSMRF 

C j 

1 L 7 ■ I 

TMP-TAUR I I > *77  3. I 60*00 
CALL  PL  ANCK  I GNII.TMP  ,BNU> 

FCL0«»0»  ( I .O.OO-M-BC  (•RNII 
27  WR|TF(A,|2) 

00  28  I-|,L 

ANG<  I ) “5  • 0*00*  OF  L 0 A T J I -I  ) 

EMUO«0CO5 • ANG  f | >*Pf/| 80 .0*00 ) 

CALL  EMNTRP ( I ) 

RRTTEC6, I3IANG( I I ,AL8UP« I ) 

01 ( I .5>«0L<  I I •FCLOMD 
E TS ( I )«8I ( t ,5) 

01 < 1 «6>»  OllliM 

00  28  J«|  |L 

PL( I I >»FMU( J)»ALBS«FP»OEXP<-<TAUA»(StGl I >*S1G< J»  )>  » 

G C I , J I ■ A ( I , J > 

AT( | ,J1-B« | ,J» 

BS< | ,JI»AC I ,J) 

AS(I«J)*A(T,J) 

PTI f . J>«8I ! , J) 

78  02  C I tJl«R( I , J> 

RFMTE  (6, 701 

WR <TE<6,3»NTF  ,TA"T  ,TAURI I ) .RNtl.ECLOUD 
RRITF(6,2»TMP0 

CALI  EMPR60(GN(I,ESHRF,FP,ALBS,TA0A,APT) 

C TAIIL*TAUT 

C CALL  OFTFCT  (CTI 

IF ( I PUNCH. EO* I 1 WRITE (3, 30) MODEL. TAUT,  APT!  l),APT<3),APT(7>iAPT<IO>. 
•APTI13) 

30  FORMAT  I | X . I | ,6X .61 F?.3 , 3X ) ) 

IF (NLYPS.EO. I IGO  TO  33 
IC»0 
C 

30  ILY«ILT*I 

IF C ILY.GT.NLYR?  IGO  TO  33 

1 C - T C ♦ I 

IF  (K.F  0.2  IRRITf  (6.701 
IF! IC.FO.I (WRITEI6.7I ) 

I F I IC.F0.2I lf-0 
TM7C«TAUR( ILY I 
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CLEM 


51 

TMP-TM7C*273. 1 60+00 

52 

CALL  PL ANCK (GNU, TMP , RNU > 

51 

EC LOUD- ( 1 •D+OO-ALBC ) »BNU 

51 

31 

CONTINUE 

IS 

no  37  |«|,L 

5i 

01  I I ,31-OLI 1 > *FCLOUD 

S7 

01  C I i 3 ) -0L ( | > •ECLOUO 

53 

Oil!  . 1 )»0I 1 I .5) 

St 

01  I 1 .21-01  I 1 .6) 

60 

IFI|LY.GT.2)ETS(!) -01(1,11) 

61 

00  32  J-I.L 

6 2 

A I 1 , J ) -AS  I | , J 1 

63 

PC  1 , J ) - A T 1 | ,J) 

63 

PI T , JI-BSI | , J) 

65 

R( I .JI-BTI 1 ,J> 

66 

ASH  ,J)-GI  | ,J) 

67 

RSI  1 . JI-GI 1 , J) 

68 

ATI  I i J ) -02  I 1 ,J) 

63 

37 

RTI T .JI-Q2I 1 ,J> 

70 

TAUU-TAIIT 

n 

A 

IF(MODFL,NE,P)GO  TO  325 

72 

TMP-0.0 

73 

GO  TO  3?6 

73 

325 

CONTINUE 

7S 

CALL  VAPOR  (HEIGHT , GNU.TMP ,99 ) 

76 

326 

TAUL-  TAUI/(CT*TMP) 

77 

WRITE (6,31 ILY  > T A UL  , TK7C , BNU , ECLOUO 

78 

WRITEI6.7ITMP 

73 

c 

AOOER  COMBINES  NEW  LAYER  WITH  OLO  RESULTANT  LAYER 

80 

CALL  AOOER I I ) 

8 1 

WR|TF(6,3)TAUT 

82 

CALL  EMPR60(GNU,ESURF,FP,ALBS,TAUA,APT) 

83 

c 

CALL  OFTFCT  I C T » 

83 

IF  I 1 PUNCH. FO, 1 (WRITE  I 3 , 30 ) MOOEL , T AUT • APT  I 1 1 • APT  1 3 1 , APT  1 7 > , APT  1 1 0 » 

85 

• APT  I | 3 I 

86 

GO  TO  30 

87 

33 

NCOUNT-NCOUNT* | 

88 

IF  I (NWAVE-NCOUNT ) .GT.O)  GO  TO  502 

83 

CALL  EXIT 

30 

STOP 

7| 

END 
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AODF  R 


1 

C 0 M P 1 L F R ( n 1 A G * 3 t 

7 

5II9R0UTTNE  A 0 0 F R (INHOMO) 

3 

IMPLICIT  RFAL»9  (A-H.O-Z) 

9 

DIMENSION  S(IA,|A,?|,T(|A,|A,2l,A(lf,,lAI,B(|A,lfc),0l(lf>,l<>t 

c 

DIHFM5T0M  AS(IA,IA),RS(IA,|6),AT(I4,I6),PT(|A,I6>,EYUPU(I6) 

A 

niMFHStON  C(1A,|A).0(|A,IAI.F(|A.IA),F(|A,|6),r,  (I6,IA),h(I6,I6) 

f 

"IMFNSIOH  Flux UP  1 I A ) ,Fl  UXON II 6 1 , ALRUP  1 I 6 ) • F YUPL ( 1 6 1 ,F YONH  < 1 6 1 

n 

DIMENSION  0|<37), E ”11(14), S|G(1AI,VM|6>,ETS(|6),E  MO  (|A) 

R 

DIMENSION  PL<IA.I6l,P(lfc,|fc),R(|6,IA>.0(|6,|A) , E T DNL l|M 

to 

FOU  I VAI.FMCF  ( 5 ( 1 > , A5  ( | ) ) , ( S ( 257  1 ,95  ( | 1 1 

1 1 

FQUI VAI  FMCF I T ( | 1 , AT ( | ) ) , ( T C 757  » ,BT  C | 1 ) 

1 7 

COMMON  S ,T  , A ,9 ,0  «QI  ,P  ,R  ,PL  ,0L  ,EMtl,FMO.»W,ETS  ,5  IG 

t 3 

common  FM,  TAIIU,FMIIO,TAIIL,TAUI, taut, alrup,  phi  ,mm,mn,l  ,ll 

I «* 

common  NT1 ,NTF ,K0 

15 

C 

* A 

00  1 1*1 ,L 

1 7 

FYUPUI l)*0|(|,|) 

IP 

FYONU ( 1 ) *Q II  ! , 7 ) 

I 9 

F Y 1 1 PI.  (11*01(1,3) 

70 

F YONL 111*0111,9) 

7 t 

II 

“> 

Ik* 

O 

C 

77 

Cl  1 ,J)*P( T ,J) 

•>7 

01  I ,J)*P( I ,J> 

■>'« 

FI  1 , JI*A5(  1 ,.)) 

?5 

FI|,J)*ATI|,J) 

7 A 

G 1 1 , J ) = 9 5 1 1 , J ) 

77 

HI  | ,.))*RTI  1 , M 

79 

T AHT*T*UL*TAUU 

79 

c 

no  99  MM  *| , HN 

30 

c 

M*MM- | 

31 

00  3 J*I,L 

37 

00  3 I *1  ,L 

33 

TA«0.'i 

39 

OO  7 K * | , L 

35 

* 

TA-TA-CI!  ,K)»E(K,J) 

3 A 

1 F 1 T .EO. J)TA*TA*I ,P 

77 

0|  1 I ,J>*  TA 

79 

** 

0 I 1 , J ) * T A 

79 

CALL  EMINV 

•ir 

00  5 J*| ,L 

•it 

DO  5 I*| ,L 

*«? 

TAO.o 

•7 

TB=n.i* 

99 

00  9 K * | , L 

•5 

* 

T A*TA*0( | ,K  ) • BIK , J) 

9A 
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00  ? 3 J ■ 1 , N 

«6 

JK  ■ JO*.) 

67 

H0L0-0M1 JK 1 

ftft 

J 1 ■ Jft*  J 

ft9 

QM( JK)a-OM( J| ) 

90 

n 

OM  C J 1 )aHOLO 

91 

29 

JaJCl  Of  ) 

9? 

IF(J-K)  20,20, ?5 

’3 

?5 

K laK-N 

99 

00  26  I at  ,N 

95 

K I aK I *N 

96 

HOLO»QM ( K I ) 

97 

JIaK l-K* J 

9ft 

OM(K 1 ) ■— OM ( J 1 1 

99 

26 

QM ( Jt ) a HOLD 

100 

60  TO  ?f» 

' ‘ " • ' 


- -r  :¥'<&*' 


SHIS  PACBf  IS  BEST  QUALITY  FJUCTI&ABLI 
rao*  oqpy  /UFsusasc  10  ddo  

E M I N V 


in| 

in* 


27  RETURN 
FND 


mis  PACE  IS  BEST  QUALITY  PRACTICAfiH 
28QI  QQPY  JTUESHISHHD  TO  DDC  ^ 


l 

fmntrp 


1 

2 

1 

9 

5 

6 
7 
A 
9 

in 
I l 
i 7 
I 7 
I 9 

»s 

I A 
I 7 
I « 
I 9 

»n 

-I 

7? 

’I 

79 

7*i 

’A 

77 

7A 

79 

in 

n 
1? 
% 7 
19 
IF 
1 A 
17 
19 
79 

• it 
*1 


compiler  ( n i a r,  * 3 i 

fhntrpimvi 

IMPLICIT  RFAI.»fl  ( A-H,n-7  I 

n IMTHS  I OH  9(|A,|A,7),TMA,IA,?),A(|6,IA),B(|A,I6).0|I|A,|6) 

n I m f m S l n m asma,ia),ps(|a,iai,at(|a,iA),bt(|a,ia),ai  p u ° (l  6 I 

0 I Mr  NS  ! 01  0 1 (17),EM|i(|A),<;ir,  (|6>.*flA>.FTS(|A),FM0<|A),FLX(9> 
n | MFIIS  ION  Pt  MA,IA),P(IA,IA),R(|A,IaI.0II6,|A) 
FQH|VAIFNCF(«;(|I.AF(|)),(S(2F7),PS(t)l 
FQI'|VAIFMCF(T(I),AT(|)I,(TI7S7),"T(|)) 

COMMON  S.T, A, 0,0.01  ,P,R,PI.  ,OL,FMU,FMn,W,FTS,S|r, 

COMMON  FM  , T All  ,FM|I0  , T AIIL  , T AU  I , T Al)T  , ALBIIP  , PH  | , MM  , M N , L » L L 
COMMON  NT  I ,NTF  ,K0 

on  r i = i ,l 

1 F I FMUO.LE  .FMIII  I ) ) c,o  TO  7 
I COMTIN'IF 

c,  0 TO  9 

7 I F I | .FO.  I)  1,0  TO  1 
I F I | . F 0 . 7 ) r, O TO  1 
IFII.ro. L ) C.o  TO  9 

| MM  * I - 7 

1 M X * I ♦ I 

CO  TO  r 
1 IMMaJ 

f M X * 9 

c,  o TO  O 

9 I M M » L -1 
1 MX  *L 

1 00  A 1 = 1,9 
A FLXI  I ) = I .OO.r-1 
00  7 ! = I MN , | MX 

| P» I ♦ | - I NN 
oo  7 J« I MN  , | MX 

IF(I.NF.J»F|XI|Pl*F|X(|PI»(EMUO-FMI|(j))/(EMUMI-FMll(J|) 

7 roNTiN'ir 
0 0 >1  J * I , N V 
ALPIIO  ( J I «0.(  A.aa 
00  R I = I MN  , I MX 
* P« I ♦ I - I «N 

« A L P 1 IP  I 1 1 »Al  p'lP  ( J I ♦FI.  X ( I P I «0  II  | , Jl 
PFTllRM 
F NO 


IBIS  PA®  IS  BfiST  QUALITY  PIUCfflCjLfiLl 
TO®  OQPY  JTJRttlsiSC  TO  flDC 


FMPR 


I 

7 

3 

6 

5 

6 
7 
6 
9 

10 
• I 
1 7 
I 3 
16 

15 

16 
I 7 
I R 
19 
70 
7 I 
’7 
73 
’6 
7« 
7 6 
77 
?R 
7« 

30 

31 
37 
'3 
99 

35 

36 
77 
3 6 
79 

•n 

tt 

97 

67 

66 

66 

•6 

67 

66 

69 

60 


C 

r 


COmpilfr  (o|ac,«j) 

SUBROUTINE  FMPR6tt(C,NII,FSNRF,FP,ALBS,TAUA,APT) 
rMPLtCIT  RFAL *R ( A-H  ,0-7 1 

I FORMAT  (/ I X .SHAMBLE  ,6X  , IOMCLOIIO  I - 0M.6X  . | OHC  l 0110  I -UP  , SX.IZHSURFA 
•F E ! -UP , 7X , | 7HSUPF Af E I -ON , 6 X , I OHTOT AL  I -ON , 7 X , | OHT 0 T At  I -UP  , ?«. 

•I3HAPPARFNT  TFMP,?X,|ZHTRANSM|SS|ON/l 
.7  FORMAT  I I X ,rs  . 7,6  I6X  ,F9.6  > I 


I 9 APP||  |976  'SECOND  VERSION* 

RE Al  TSS , TST , Y ,03 

0 I MENS  I ON  S(l6,|6,?>,T(16,|6,2>,A(|6,l6),B(|A,l6>,(}|(|6,l6t 

DIMENSION  AS(I6,I6),BS(|6i|6I,AT(I6,I6),B.T(|6,I6),ALBUR(I6) 

n I Mr  NS  I ON  01  I 3 7 I • E MU  I I 6 I , S I C,  (|6>.AMI6>,ETS(|6I,EM0(|6),Y(|6> 

d|MFNS|ON  PL<l6,|6t,PII6, 161, R(|6,I6>, 0(16,161, APTII6I 

niMFMSION  0(16,161,  07 ( | 6 , I 6 I ,EYUP ( I 6 ) ,FT0N ( | 6 I , TAUR (SO  I 

F0I!|VAIF*ICF(S(|J,AS(|II,(S(2S7I,BSII)I 

rO"IVAI  ENCF ( T ( | I , A T ( | ) I , (T(7S7I  ,BT(  | 1 I 

COMMON  S , T , A ,R  ,0 ,0  I ,P  ,R  ,P|.  ,01  , E Mil  , F M 0 , * , E T S , S 1 C, 

rOMMON  fm,  TAUII.EMUO.TAUL  , T AU I , TAUT . »L6UP  , PH | , MM , NN , L • LL 

COMMON  NT | , N T r , K O 

OO  7 | » | ,L 

7 E TUP ( I 1«SIB(  I I • F MU ( | t»FSURF»OfXP(— TAHA»51G(  I • » 

IE ( Al  BS.ro. 0.0*00  I 00  TO  IP 
00  S I = I ,L 
00  S J« | ,L 
TA  ■ P.O+Oll 
00  6 K«t  ,L 

6 T A * T A — P|.  ( I , K I • B S ( K , J » 

IF(|.EO.J)TA»l.0*C0*TA 

S oil  ,.I|,TA 
CALL  E«|NV 
00  7 |«| ,L 
00  7 J*|,L 
T A *0 • 0*00 


TB»0»0*00 

00  6 K ■ | , L 

TA»TA»«T  ( I , K I • O ( X , .)  ) 

6 TB«TB*bs I I ,K I *0  I E , Jl 
B(  | , J I * T A 

7 PI  I ,JI«TB 

00  9 | « | , L 

00  9 J«  | ,L 
TA  ■ 0.0*00 
00  6 K • | , L 

6 TA*TA*n(l,K)*R| ( K , J I 
9 R I I , .1  | ■ T A 

10  I E ( Al  RS  .NE  .1. 0*1:0  » So  TO 
00  II  l«l ,L 
00  Cl  J 6 1,1. 

P I I , J ) -RSI  | , .1 1 


I 7 


i 


\ 

\ 


1 
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SHIS  Pi.dE  IS  BEST  QUALITY  PRACTICABLE 
IBM  OQPY  /UF5HLSHM)  TO  DOC  - 


i > 

fMPP 


¥ 


r-i 
E ? 

93 

r*t 

•■5 
r a 
'7 
r 4 
t Q 

ao 

4? 

43 

44 

44 
47 
4* 
49 
70 
7| 
7? 
* 3 

79 
79 
74 
77 
74 
79 
40 
"I 
4? 
"3 
49 
95 
44 
47 
9 4 
49 
90 
9 I 
97 

93 

94 

95 
94 


PIT,  I 1 * 4 T ( f , J > 

P ( I , J I *0 .D*jP 

II  pl  ( i • j ) *p.n ♦«'!' 

| 7 no  | 4 I » | ,1 

TA»n.n*ni 
TR=o.o*ro 
no  7 3 I*  I ,L 

r»*sir,  (i  i i 

t a = t»  *f a *n ( i » j i • f v up ( j t 

I 7 TB  * TP  ♦ FA*P(I,J)*Q|(J,B) 

0 | ( ! ,7  I »TA 

I 4 0 | (|  . I I > *TB  *0  I (|,4I 
no  14  I * I , L 
no  14  i«i  ,l 

T A *f.  • 0 *00 
00  15  F»  | ,| 

| 5 T A * T A ♦ P ( I , J I • P L I K • J I 
14  P(  I ,.I)*TA 
00  14  | « | ,L 

TA»0.0*m 
T 4 «0 • 0 ♦OO 
no  1 7 .1*1 . u 

F4»516IT»»«*IJ)»FMU(J|/«Mt) 

ta*ta*fa*p( i , j i *r vup i j i 
I 7 TB*T4*FA*R(  | , J I • 0 | IJ.5I 
Oil T , IP7-TA 
Q I I I ,4  1 » TP ♦ 0 I ( | ,5  I 
0| ( | ,9  1 *0| I | ,4 ) *01  I | , 101 
OII|.|4)«IO|lt,m-IPH|»OlHtn/fTS(|l 
FT5 I I I >0 II  I , I I > 
l«  01  I I , I 7) »a|  ( I ,7  I *01(1,11) 
no  19  T * I , I 

0 I ( I . I 7 ) *0 | ( I ,7  I /E5IIRF 

1 9 oi(i,i*n»((i . 990*00* gnu  )/(mo6(i.n*oo*((i.  1 9340-05  *r,Mu**  3.  o*oo ) 
•/I09l|,|7)|))l 1-773.140*00 
WR I TF I 4 , I ) 

PAn5«i.n* or /(fp*7. 70*07) 
no  70  1*1,14 

ANO  *5.00*r J*OFLOAT ( l-l  I 
FMIIO"OF05<AHG*PAOSI 
CALL  EMNTRPII4) 

APT  ( t ) *AI  Blip  ( I 9 ) 

?r  PR  TTF  I 4 , 7 I ANfi  , At  4UP  ( 4 ) , Al  BlIP  I I I ) . »LBllP  ( 7 ) , AL41IP  ( 10)  , AL9|)P(  9 I , 

• A L 4 ' IP  ( 17)  ,ALBIIPI  14)  ,ALBIIP(  14) 

RETURN 

FNO 
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Ibis  pase  is  bbst  quality  nufiCffillU 

FROM  COPY  FURBISHED  10  DDC  


I n 


SlIXRrtllT  t HF  in 
niMFf'sinn  nuMMvoo) 

PRINT  I »»n 

|*TP  FORMAT  ( I H | , • • • • T H F F0L1.0»I1MG  IS  THF  PROPFP  INPUT  DATA* 
♦'  TO  ftn  *ITH  THE  SURSF01IFNT  OUTPUT', /I 
no  | l = | » 3 R 

FFAD  ( S , 1 0 I , FNO  = ? I ( n IIMM  v ( ,|  ) , j » | , ?n  | 

»R!TF  (6.IPI)  ( DUMMY  ( K I , F » | , ?f# ) 

I rOMTIMUF 

? rOMTIMUF 

PRINT  in? 

10?  rORI'AT  ( I H ? , ' FMO  INPUT') 

I »■  I FORMAT  I ?n  A H ) 

RFTIIPI4 
F NO 


THIS  PACE  IS  BEST  QUALITY  PRACTICABjyi 
raOli  COPY  TlfSHiSHH)  IODDC  


PLANCK 


roMPtLfR  inuciji 

SlIHROllT  I >JF  PI  «MCK  ( GMII  , TMP  , null  ) 
r»  0 1 1 p L E "RECIMON  PNM,TMP,RNII,0FXP 
iriT«p.i.T.x1n*i>i  ir,n  rn  i 

/ ( r>F  X P ( 1 ,<mn*OP»fiNli/TMP)-|  . p ♦ 00  I 

r, 0 TP  ■> 

I p nm  * r • n *r>r 
» RETURN 
Fim 


MB  IS  BBT  waznntm&m 

JS»PY?v3ra»lsHB)TOl^C  — 


VAPnp 


I 

7 

7 

9 

9 

A 

T 

P 

9 

1" 
I I 

I 7 
t 1 
19 
|f 
I A 
t 7 
I « 
I 9 

■»n 
7 1 

77 

71 

79 

♦ C 

76 

77 
79 
’9 
IP 
71 
'? 
n 
79 

IS 
7 4 
77 
79 
79 

90 
9 I 
9? 

91 
•9 
95 
94 
97 
9 9 
99 
«0 


CUHPiLro  < n 1 a s * 3 > 

SMBpniiT  r HE  vaporIhT.gniio.tan, MOOFL) 

DOUBLE  PRECISION  TAII,  GNI>n  ,HT 

C FOP  r Al  f ML  A T ! NG  OPTICAL  "rpTH  (TAMA  I OF  FATEP-VAPOR  (D"E  TO 

C C0M7INMIIM  ABSORPTION  1 M TMF  Fort  TO  |7S0  CM-|  RANGE)  PFTwFFM  THF 

C I1  A Sr  or  a noun  (HFIGHT*0|STAIICF  1*1  IN  9f  T 4F  rM  BASF  ANO  SURF  A C F ) 

C AN"  Thp  BUPFACF. 

C MOOFI  I IS  TMF  AFcpl  TPOPICAI.  ATmosPhFP  I C POOFL 

C M 0 O F I 7 (9  TMF  Ar  CPL  MIOLATITMPF  9UMMFR  POOF  L 

C MOOFI  1 IS  TMF  A F C R L M!0|  ATITMOF  w-INTFP  moDFI 

C MOOFI  9 (9  THF  1 9 A 7 IIS  STAMIjapo  ATMOSPHFPF  MOOFL 

C ••  Sr E rFL9Y (ROHFRT9  ♦ 9 I 0 F R M A M ( | 9 7 A ) FOR  A9S0RP.  MODFL  •• 

C 

C MOOFI  =99  A S 9 m M F S 9ATMRATI0N  AT  TFpp* 1 T AU  * 

0 (mf  MS  I nil  7l|i’),»VR(|n),TP(ir),WRM(|0,9),TRH{|r,9) 

OATA  //'•.,!. ,7. ,3.  ,9. ,9., A. ,7. ,9. ,9./ 

OATA  APH/|9.,|3.,9.3,9.7,7.7.I.S.»99,.97,.7S..|?i 

• I9.,9.3,9.9,1.1,|.9,f.,.6|,.37,.7l,.l2, 

• 3*9, 7.9, I .9, |. 7, .44, *39, .21, .099, .03^, .016, 

• 9. 9,9. 7, 7. 9, t. 9, I. 1, .69,. 39 17, .094/ 

OATA  T PM / 3f l • ,799.  , 7 9 ft . ,799.  ,277.  ,770. , 7 A9 . ,757. ,7S0.  ,799.  , 

• 799.,79».,?PS.,779.,773.,?A7.»76*.,?SS.,79B.,292., 

• 777. ? , 7 6 9.7 , 7*9. 7 , ?4 I . 7 , 70S. 7 ,799.7 , 793.7 , 737.7 , 73  I . 7 , 2 29.7 , 

• 709.|,7A|,a,?7S.I  ,749.7,767.7, 799.7 ,7<(9. 7, 792. 7,736. 2, 729. 7/ 

I rOPMMI/1  ••••  THF  VALMF  GNM  = ',F|O.A,'  IS  OUTSIDE  THF  DATA  RANG 

• F F 0 9 CONTINUUM  ABSORPTION  DATA.  TAIJA  SFT  * 0«0  •••/) 

9 FORMAT!/'  •••  ri OMO-RASF  HFIGHT  of  ',E|7.9,'  IN.  EXCFFOS  INTERNAl 

• OATA  RANGE  ...HFIGMT  proHCFO  TO  9.0  F M ♦ ♦•/.) 

7 FOPMATI/IX,'  ATMOSPHTRIC  MOOFL  tor  9IIR-CI OMO  ABSORPTION  » ',12,'. 

• OPTICAI  OF  n T M ( TAMA  ) = •,FI2.S,'  CLOUO  RASE  HFIGHT  » * , E I 2 ♦ S , • 

• KM'/) 

7 FORMAT!/'  ♦♦♦  MOOFL  NO.  ',13,'  IS  UNOFFINFD.  TAMA  SFT  TO  ZERO*/ 

• ) 

G N M ■ S N G | ( GNIlO  I 

1 F IGUM.I.t  .900.0. OR  .GNU.GT.  I 2SP.  )G0  TO  3 
GO  TO  9 

1 FR  ?TF I A , I I GNU 
GO  TO  A 

« HE  I 6HT*9NGt  ( HT ) 

IF  IMP  | GHT  .|.F  .9.  >G0  TO  S 
FR t TF | 6 , 9 » HFIGHT 
HE  t GHT  »9 .0 

5 CSO«| .79F-72*  C | . ♦ 1 1 36 ♦ »F X P I -7 . 9 7F -3 »GNU I ) 

C A YC*P .0 

IF(H0DFI.,E0. 99)60  TO  70 

IF  I MOOFL  .C,F  . I . ANO.MOOFL  «LF  . 9 ) GO  TO  9 

6 CAYC«lUn 

IF IMOOFI ♦EO.PIGO  TO  99 
WRtTF(6,7)M0PFL 
GO  TO  99 


IBIS  PAQS  XS  BEST  QUALITY  FRACTICABLS 
IBM  OQEX  JMWISHH)  TQVO&  — ^ 


VAPOR 


r I 

c 1 

F3 
r9 
tt 
CA 
* 1 
r A 
c 9 
AO 
A I 
A 7 
A 3 
A A 
AF 
A A 
A7 
AB 
A9 
7r 
T I 
7» 
73 

7*1 
7f 
7 A 
77 

7 A 
79 
AO 
«l 
«? 
»3 
A«l 
AF 
AA 
A 7 
Afl 
«9 


p no  1 < t « i , ic 

Y'R  I I t *#BH  ( ! , MODE  L ) 

I r*  T P t t ) ■ T i»  H ( | , Monf  L 1 
tjM»s|  ♦ IFtxihfifhT) 

Tl  AFTM  F I X t (HFIfiHT-FLOATlFMX-l  ) Min.  ) 

IF(IHE»r,MT-(FLnAT«HMX-ll*.|»FLnAT(llAST)>l.r,F.f».'1S>IIAl;T*ILAST*| 
I M X « I O 

CAvr*u*o 

on  Ail  M*t  ,MMX 
! F t N.FO.NMX  ) I MX  = It.  AFT 
|FM«X.FO.rir,n  TO  99 
0 7*7 ( M ♦ 1 )-7(f) 

A | * < WR  t N I *7  ( N*  t I -BP  (*(♦  ( I • 7 t H I I / n 7 
A 2«  ( HR  » M*  | |-»IA(N)  I/07 
A|*ITRIH)»7IN*|  I “TP  t N* I I » 7 I N > ) / 0 7 
P ?* I TR ( H* I ) -T9  < N ) ) / n 7 

nn  fa.  1 = 1 . 1 m x 

7C«FLO*T(M-| )♦( ( .t»FL0ATI | I ) - .OF ) 

»*A|«A7*7r 
T * A ! *B7»7r 
PIM9.5AF-SM** 

» H »fl  * A « 3 . 3 A F I A 

FACT*C*n«PtM"M*c»FXPI |»or.»t I |./T)  -.003370)1  • I.F5 
f A Y C * C A Y C ♦ F AFT 
= F rni'TINIIF 
A*»  rOIITINMF 
99  TAUA*CAVf». I 

BRITF{A,7|MonFI.  ,TAHA,HF|r,HT 
TAIImOBI  F ( T All  A ) 

FO  TO  ion 
70  T=FNFL(TAU' 

FF»A.|°A  •FXPtF.P%3FOM3.AA3E-3-tl.7T))) 

w«?1 A.F.FS/T 

PIMA  .FAF-F»T»*  ••! 

WH?0***3.3AEI  A 

r A YC*CFn»Pw«WH?o*FXP t I BOP. • < I I . 7T  ) -.0O337B1)  • |.EF 
TAIMOBI.  F 1 CA  YF  ) 

lijr  return 
f Nn 


B B8ST  9UALrTy  P3W,CTIG4BL1 
COPY  7URSUSHBC  roicc  vw“* 


SAMPLE  INPUT  DATA 

**  THE  kolLO^ING  IS  THE  PKOPER  INPUT  DATA  TO  GO  WITH  THE  SUBSEQUENT  OUTPUT 

LO  O.IOOUOO  0.000000  30.000000 

1 2.0  +00.00001 

0.2500  -0.7500  -1.2500  -1.7500  -2.2500  -2.7500  -3.2500  -3.7500  -4.2500  -4.75o0 

5.2500  -5.7500  -6.2500  -6.7500  -7.2500  -7.7500  -6.2500  -8.7500  -9.2500  -9.75oO 
1 

32  .1UOO00+004  .115078+003 
.67070903723256+000  0 

.10514996324114+001  1 

« 2657 54761b07bl +001  2 

.35790199403719+001  3 

,4Oo0O135473 38 1+001  4 

.43356235234477+001  5 

,44/1632204*415+001  6 

. 4444 J36462u930 +001  7 

.429t»7e>2l537b35  + 001  0 

,4OoO02465oo076+OOl  9 

.37314599004131+001  10 

,334457O40«2b67+OOl  11 

.30073792306375+001  12 

•2ou904 27976294 +001  13 

.22149295911243+001  14 

. l63692o5002  706  + 001  15 

.I4rt409o5736291+001  16 

.11647350560741+001  17 

.60304740775969+000  10  ^ 

.64104702353509+000  19 

,4420664b072O24+OOO  20 

,2047o653Ole6O9+OOO  21 

.16727432440854+000  22 

.07290112986402-001  23 

.40440299357566-001  24 

.16154620850324-001  25 

.57360106061632-002  26 

.17932424533112-002  27 

,499369b0377115-003  28 

.12508217919023-003  29 

,204l74254lo692-OO4  30 

.57262299763976-005  31 

end  Input 
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SAMPLE  OUTPUT  LISTING 


Mis  PAGE  IS  BEST  QUALITY  FRACTICABUI 
00 PY  EURMISiffiD  XO  dqq 


M O m X 
£ 0»<(M 
* »o  ^ p- 
X X X n 
trt  CV  0 O 

*o  a m o 
4 <0  p-  o 

4 m X O 
4 CC  O O 


4 *+  4 p* 
o O'  rsj  « 
'M  41  m 4 
\'  O'  * X 
41  <V  N CM 
X O'  "\J  O 
^4  4 p-  o 
p -»  x o 
» <\i  H O 


p*  4!  o «-• 

-«  cm  *o 

!*>  x o cm 

'-1  7>  O tJ 

J\  K if*  o 

■O  C\J  X 41 
*M  4 O CM 
X X1  P-  — 
41  7 4 C 
41  o 4)  o 
in  fw  c 


p ~4  0s  O 
4 o o ■** 
41  cm  4 »o 
•O  O'  O O' 
*4  N OtT1 
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